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1.0  Executive  Summary 

The  U.S.  Air  Force  is  interested  in  advanced  micropropulsion  systems  for 
nanosatellite  applications.  One  of  the  challenging  areas  in  the  development  of 
micropropulsion  systems  is  the  design  of  micronozzles.  Such  micronozzles  are  needed  to 
generate  very  small  thrusts  to  enable  microsatellites  to  counteract  the  disturbances 
produced  by  solar  pressure,  drag,  or  the  distortion  in  the  gravity  field  as  well  as  to 
reposition  microsatellite  instruments  according  to  mission  requirements.  Thus,  both 
experimental  studies  and  modeling  are  needed  for  the  successful  design  of  micronozzles. 

During  the  course  of  this  research  effort  gas  flows  in  microchannels  and 
micronozzles  were  studied  both  experimentally  and  numerically.  For  the  experimental 
study  a  flow  visualization  system  was  built  and  used  to  study  gas  flows  in  microscale. 
Gas  velocity  measurements  in  microscale  were  conducted  using  both  Laser  Induced 
Fluorescence  technique  (LIF)  in  conjunction  with  Image  Correlation  Velocimetry  (ICV) 
and  Molecular  Tagging  Velocimetry  (MTV)  technique.  Maximum  deviation  of 
experimentally  measured  average  velocity  via  LIF-ICV  technique  from  estimation  based 
on  the  gas  mass  flow  rate  was  12.4%. 

Three  different  approaches  were  utilized  in  the  numerical  study.  Continuum 
computational  fluid  dynamics  was  used  to  study  gas  flows  in  microchannels  and 
micronozzles.  For  flows  in  micronozzles,  effects  of  geometrical  scaling  down  and 
different  gas  propellants  were  studied.  A  correlation  for  specific  impulse  involving  throat 
diameter  and  throat  Reynolds  number  was  developed.  For  flows  in  microchannels, 
utilization  of  slip  versus  no-silp  boundary  condition,  compressibility  and  rarefaction,  and 
aspect  ratio  effects  were  studied.  An  analytical  correlation  was  developed  to  predict  the 
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normalized  friction  factor  variation  with  change  in  pressure  ratio  for  both  incompressible 
and  compressible  slip  flows.  For  compressible  flow  in  a  microchannel,  assumption  of  the 
flow  to  be  two  dimensional  was  found  to  result  in  prediction  of  15%  -  45%  higher 
velocities  and  7%  -12%  lower  normalized  friction  factor,  then  the  corresponding  three 
dimensional  solution. 

Direct  simulation  Monte  Carlo  (DSMC)  method  was  used  to  study  low  Reynolds 
number  flows  in  a  conical  micronozzle.  The  nozzle  exit  thrust  simulated  by  DSMC 
showed  good  agreement  with  experimental  data  from  open  literature  at  throat  Reynolds 
number  larger  than  10,  where  the  flow  is  in  slip  regime.  For  throat  Reynolds  numbers 
below  10,  the  DSMC  predicted  larger  thrust  than  the  experimental  data. -The  DSMC  was 
also  used  to  study  propellant  gas  temperature  effect  on  the  generated  thrust. 

A  Unified  Flow  Solver  that  utilizes  hybrid  approach  using  deterministic 
Boltzmann  solver  for  highly  non-equilibrium  flows  at  high  Knudsen  number  and 
continuum  (Euler  or  Navier-Stokes)  solvers  for  low  Knudsen  numbers  was  tested  and 
demonstrated  for  gas  flows  in  microscale.  Test  cases  included  both  microcnannel  and 
micronozzle  flows  in  continuum,  slip,  transition,  and  free  molecular  flow  regimes. 

Six  students  were  involved  in  this  research  effort.  The  outcome  of  this 
involvement  was  the  production  of  Masters  thesis  by  V.  V.  Gadepalli  entitled  “Numerical 
Simulation  Of  Gas  Flows  In  A  De-Laval  Micro  Nozzle.”  Production  of  a  Ph.D.  thesis  by 
S.  Gokaltun  entitled  “Particle-based  Simulations  for  Rarefied  Gas  Row  in 
Microgeometries”  is  in  progress.  The  expected  outcomes  of  this  thesis  are  adaptation  of 
Lattice  Boltzmann  Method  (LBM)  for  gas  flows  in  microscale  and  development  hybrid 
solver  for  modeling  microscale  gas  flows  that  combines  DSMC  and  LBM  approaches. 
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Major  conclusions  from  the  thesis  by  V.  V.  Gadepalli  are  as  follows:  1.  Employed 
continuum  model  was  able  to  simulate  nozzle  exit  thrust  which  is  in  good  agreement  with 
the  experimental  data  except  at  very  low  Reynolds  number;  2.  The  model  predicted  closer 
values  to  the  experimental  data  even  for  the  flows  in  transitional  flow  regime  at  the 
nozzle  exit.  This  indicates  the  model  can  also  be  used  for  transitional  flows;  3.  A 
correlation  for  specific  impulse  involving  throat  diameter  and  throat  Reynolds  number 

was  developed  as  /  =  29.54  •  Re3  D~2M .  The  correlation  can  predict  Isp  with  a  marginal 

error  up  to  10  %  of  the  actual  value;  4.  In  a  study  of  gas  propellants  is  was  found  that 
helium  proved  to  be  best  among  nitrogen,  argon  and  carbon  dioxide  for  throat  Reynolds 
number  varying  from  5  to  2000  for  the  three  throat  diameters  investigated.  Nitrogen 
followed  helium.  The  performance  of  argon  and  carbon  dioxide  varied  depending  on  the 
throat  diameter  and  the  throat  Reynolds  number. 
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2.0  Introduction 

The  U.S.  Air  Force  is  interested  in  advanced  micropropulsion  systems  for 
nanosatellite  applications.  One  of  the  challenging  areas  in  the  development  of 
micropropulsion  systems  is  the  design  of  micronozzles.  Such  micronozzles  are  needed  to 
generate  very  small  thrusts  to  enable  microsatellites  to  counteract  the  disturbances 
produced  by  solar  pressure,  drag,  or  the  distortion  in  the  gravity  field  as  well  as  to 
reposition  microsatellite  instruments  according  to  mission  requirements.  Thus, 
micronozzles  should  be  capable  of  producing  highly  accurate  impulses  or  steady-state 
thrust.  Consequently,  the  design  of  micronozzles  is  very  important.  It  has  been  shown 
that,  in  many  instances,  conventional  fluid  dynamics  approaches  are  not  valid  in 
microscale.  The  fundamental  flow  processes  in  microscale  are  currently  not  well 
understood.  Thus,  both  experimental  studies  and  modeling  are  needed  for  the  successful 
design  of  micronozzles. 

The  main  objectives  of  this  research  effort  were  to  develop  and  demonstrate  a 
flow  visualization  technique  capable  of  measuring  gas  flow  velocities  inside  a 
micronozzle,  to  characterize  flow  inside  a  microchannel  and  several  converging  and 
diverging  micronozzles,  and  to  compare  experimental  flow  measurements  with  the 
predictions  of  existing  microscale  flow  models  with  the  aim  of  validating  and  improving 
these  models. 

For  experimental  measurements,  application  of  Image  Correlation  Velocimetry 
(ICV)  technique  to  velocity  measurements  in  microscale  gas  flow  has  been  extended. 
Image  correlation  velocimetry  is  a  spatial  image  correlation  technique  allowing 
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estimation  of  the  displacement  (and  consequently  velocity)  vectors.  This  technique  has 
been  successfully  used  by  Genderich  and  Koochesfahani  (1999)  and  Fielding  et  al. 
(2001)  among  others,  for  velocity  measurements  in  conventional  scale. 

Recently,  Ayon  et  al.  (2001)  successfully  demonstrated  the  deep  reactive  ion 
etching  (DRIE)  technique  for  fabrication  of  silicon  supersonic  micronozzles.  They  also 
tested  the  thrust  performance  of  several  nozzles.  However,  no  detailed  measurements  of 
flow  field  inside  or  outside  the  nozzles  were  performed.  For  velocity  profile 
measurements  of  gas  jets  produced  by  micronozzles,  several  velocity  measurement 
techniques  have  been  used,  such  as  laser  doppler  anemometry  (LDA),  particle  image 
velocimetry  (PIV),  molecular  tagging  velocimetry  (MTV)  (Lempert  et  al.,  2002),  and 
planar  laser-induced  fluorescence  (PLEF)  (Kato  et  al.,  1991). 

For  velocity  profile  measurements  inside  a  channel,  Meinhart  et  al.  (1999)  applied 
the  microscale  particle  image  velocimetry  (micro-PIV)  technique.  The  same  technique 
was  applied  by  Wereley  et  al.  (2002)  to  study  flows  inside  a  micronozzle.  However,  they 
used  liquid  to  operate  the  micronozzle  because  of  difficulties  in  seeding  the  gas  flow. 
Because  of  the  unavailability  of  suitable  seeding,  no  PIV  measurements  in  gas  flow  in  a 
microchannel  have  been  reported  to  date.  MTV  has  also  been  used  to  measure  velocity 
profiles  in  liquid  flows  in  microchannels  (Paul  et  al.,  1998;  Maynes  and  Webb,  2002; 
Thompson  et  al.,  2002).  To  the  best  of  the  authors’  knowledge,  the  application  of  ICV  to 
gas  flow  in  microchannels  has  not  been  attempted.  This  technique  has  the  advantage  of 
not  requiring  seeding  particles  to  operate,  and  have  been  demonstrated  for  microscale 
flows  applications  in  this  report. 
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From  the  modeling  standpoint,  there  have  been  a  number  of  efforts  to  model  fluid 
flow  in  microscale.  Some  of  the  recent  relevant  publications  include  Koplik  et  al.,  1995; 
Gad-el-Hak,  1998;  Beskok  and  Kamiadakis,  1999;  and  Simek  and  Hadjiconstantinou, 
2001.  In  summary,  microscale  flow  regimes  can  be  characterized  by  the  Knudsen  number 
(Kn): 

•  Kn<  0.001,  continuum  flow  (the  continuum  model  can  be  used  reasonably); 

•  0.001  <  Kn  <  0.1,  slip  flow  (the  continuum  model  with  slip  boundary  conditions 
can  be  used); 

•  0.1  <  Kn  <  10,  transition  flow  (Direct  Simulation  Monte  Carlo  (DSMC)  method 
can  be  used);  and 

•  Kn  >  10,  free-molecule  flow  (DSMC  method  can  be  used). 

In  this  report  continuum,  DSMC,  and  hybrid  solvers  were  used  to  check  their 

applicability  to  various  flow  regimes  and  perform  parametric  studies  of  gas  flow  in 
microscale  to  understand  the  physics  of  such  flows. 

3.0  Research  Results 

3.1  Experimental  Study 

3.1.1  Experimental  Set-up 

A  schematic  diagram  of  a  test  set-up  for  studying  gas  flows  in  microchannels  and 
micronozzles  used  in  this  research  effort  is  shown  in  Figure  1.  The  set-up  consists  of  a 
gas  supply  cylinder,  a  mass  flow  controller,  an  acetone  bubbler,  a  test  cell,  manometers,  a 
dump  tank,  and  a  vacuum  pump.  The  set-up  operates  as  follows.  The  test  cell  is 
evacuated  by  the  vacuum  pump,  which  is  connected  to  the  cell  trough  the  dump  tank/cold 
trap,  the  valve  between  the  pump  and  the  dump  tank  allows  pumping  rate  regulation.  The 
dump  tank/cold  trap  is  used  to  capture  all  acetone  vapor  from  the  working  gas  before  it 
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enters  the  vacuum  pump.  The  working  gas  is  supplied  from  the  gas  cylinder,  its  flow  rate 
is  controlled  by  the  mass  flow  controller.  Seeding  of  the  working  gas  with  acetone  (for 
LIF  detection)  is  done  in  the  acetone  bubbler.  The  bubbler  is  placed  into  the  water  bath  to 
allow  temperature  control  and  thus  acetone  concentration  control.  The  bubbler  can  be 
bypassed,  thus  acetone  concentration  gradients  could  be  introduced  in  the  working  gas. 
The  bypass  valve  is  the  3-way  solenoid  valve,  with  flow  switching  frequency  controlled 
by  a  data  acquisition  system.  Two  absolute  pressure  transducers  are  monitoring  pressure 
at  the  test  section  inlet  (microchannel/micronozzle)  and  inside  the  test  cell.  One 
thermocouple  is  installed  into  the  acetone  bubbler  and  the  other  into  the  dump  tank/cold 
trap  for  temperature  monitoring  of  these  devices.  All  test  set-up  instrumentation  is 
connected  to  a  data  acquisition  system,  which  facilitates  control  of  the  set-up  operation 
and  data  collection  via  LabView  code  developed  in  house. 

Mass  Flow  Manometer  1 


Figure  1.  Schematic  diagram  of  the  experimental  set-up. 
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The  test  cell  is  shown  schematically  on  Figure  2.  The  test  cell  consists  of  the  cell  and  the 
cover,  which  were  machined  out  of  aluminum  alloy.  The  cover  has  three  ports  for  gas 
inlet,  pressure  monitoring,  and  vacuum  line.  The  cover  is  sealed  onto  the  cell  with  the 
help  of  six  bolts  and  an  o-ring.  The  gas  inlet  port  is  connected  to  a  cylindrical  plenum  that 
provides  uniform  flow  condition  at  the  microchannel  inlet.  For  optical  access  the  cell  has 
three  quartz  windows,  one  on  the  bottom  and  two  on  the  sides  as  shown  on  the 
photograph  in  Figure  3. 


To  the  Pressure  To  the  Vacuum 
Gas  Inlet  Transducer  System 


Cover 


1  t  t 

VD _ □ _ H_ 

bn — i — r— - cd 


1 _ I _ u 

■ -  Plenum 

MicroChannel 

Figure  2.  Schematic  diagram  of  the  test  cell. 
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Figure  3.  Test  cell. 

Test  apparatus  validation  was  performed  by  conducting  flow  tests  with  pure 
nitrogen.  In  these  tests  the  pumping  rate  was  adjusted  with  the  valve  between  the  dump 
tank  and  the  vacuum  pump,  and  then  the  gas  flow  rate  was  changed  with  the  help  of  mass 
flow  controller  and  pressures  at  the  channel  inlet  and  outlet  were  measured  with 
manometers  1  and  2.  Quartz  microchannel  with  square  cross  section  was  used  in  the  tests. 
The  inside  dimension  of  the  channel  is  500  pm,  and  the  length  of  the  channel  is  12.7  mm. 
Reynolds  number  for  conducted  tests  was  in  the  range  of  30  to  300,  while  Knudesen 
number  was  between  lxlO'3  to  8x1 0‘3. 

Measured  pressure  drops  across  the  channel  for  different  gas  flow  rates  were 


compared  with  correlations  for  laminar  flow  in  rectangular  ducts  (Shah  and  London 
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1978,  Shah  and  Bhatti  1987).  This  comparison  is  possible  because  the  flow  in  the  channel 
is  in  continuum  regime  as  indicated  by  low  Kn  number.  Following  correlations  were 
used.  For  fully  developed  flow  in  a  square  duct: 

/Re  =  14.227  (1) 


where  /  = 


is  the  Fanning  friction  factor,  and  Re  is  the  Reynolds  number. 


The  non-dimensional  hydrodynamic  entrance  length  for  the  square  duct  is 
+  Lhy/Dh 

Lh  =  — — - =  0.090 ,  where  Dh  is  the  hydraulic  diameter  (equal  to  duct  height  or  width 

Re 


for  the  square  duct).  In  the  entrance  section  of  the  square  duct: 

/wRe  =  18 

The  pressure  drop  from  duct  entrance  to  exit  is: 


(2) 


Ap=^(/„Re-Lt,  +  /Re(L-zJ)  (3) 

Uh 

where  L  is  the  length  of  the  duct. 

The  comparison  between  the  measurements  in  the  square  microchannel  and 
square  duct  correlations  is  shown  in  Figure  4.  Two  duct  hydraulic  diameters  of  500  pm 
and  550  pm  were  used  in  the  correlations.  Figure  4  shows  that  experimental  data  falls 
between  the  curves  predicted  by  correlations.  The  discrepancy  could  be  attributed  to  the 
fact  that  channel  cross  section  has  rounded  comers  as  revealed  by  studying  the  channel 
under  a  microscope,  while  correlations  were  developed  for  a  duct  with  perfect  right 
angles.  Also  nonuniformity  of  the  channel  cross  section  along  the  channel  length  could 
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affect  apparent  hydraulic  diameter,  the  manufacturer  specifies  channel  wall  tolerance  of 
±150  pm  for  this  channel  with  nominal  wall  thickness  of  250  pm. 


Figure  4.  Comparison  of  experimental  data  with  square  duct  correlations. 

3,1.2  Micro  Gas  LIF  System 

Existing  micro  PIV  system  was  upgraded  and  modified  to  enable  this  system  to 
acquire  LIF  images  in  acetone  seeded  gaseous  flows.  The  Nd:YAG  Solo-120  PIV  laser 
was  upgraded  to  have  ultraviolet  (266  nm)  output  to  allow  acetone  excitation.  Proper 
filters  were  installed  in  the  Nikon  Eclipse  TE2000  microscope  to  allow  acetone 
fluorescence  observation.  An  image  intensifier  and  a  synchronization  unit  were  added  to 
the  system  to  allow  LIF  image  acquisition  by  the  CCD  camera.  LIF  add-on  software 
module  was  installed  on  the  PC  controlling  the  system  to  allow  LIF  image  processing. 
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The  diffraction  limits  of  the  microscope  with  lOx  and  40x  objectives  were 

X 

calculated  as  (d=2.44-M - ,  Meinhart  et  al.  1999)  12.9  pm  and  21.6  urn 

2-NA 

respectively,  which  corresponds  to  the  point  size  projected  onto  CCD  sensor.  These  point 
sizes  are  easily  resolvable  by  our  CCD  camera  with  pixel  size  of  6.45  pm  and  pitch  size 
of  6.45  pm.  When  these  values  are  projected  back  into  the  flow  they  correspond  to  spatial 
resolution  of  1.29  pm  and  540  nm  for  lOx  and  40x  objectives  respectively.  When  pixel 
binning  is  used  the  spatial  resolution  is  limited  by  the  CCD  sensor  pixel  size.  For  2  by  2 
binning  the  spatial  resolutions  for  lOx  and  40x  objectives  are  2.58  pm  and  645  nm 
respectively.  For  4  by  4  binning  the  spatial  resolutions  are  5.16  pm  and  1 .29  pm  again  for 
lOx  and  40x  respectively. 

3.1.3  Image  Correlation  Velocimetry 

Image  correlation  velocimetry  (ICV)  procedure  was  developed  using  MatLab. 
The  procedure  uses  algorithm  of  Gendrich  and  Koochesfahani  (1996).  The  developed 
MatLab  script  will  be  used  to  analyze  LIF  images  acquired  by  micro  gas  LIF  system. 
Dantec’s  software,  called  Flow  Manager,  for  gas  LIF  system  acquisition  and  control, 
provides  MatLab  link  so  that  custom  image  analysis  procedures  can  be  easily  performed 
using  MatLab  scripts.  The  results  of  such  custom  image  analysis  are  conveniently  stored 
in  the  same  database  with  originally  acquired  images  simplifying  analysis  process  and 
reducing  possibility  of  human  errors  while  working  with  large  number  of  images. 

Image  correlation  velocimetry  is  a  spatial  image  correlation  technique  allowing 
estimation  of  the  displacement  (and  consequently  velocity)  vectors.  The  fundamental 
assumptions  of  this  technique  are  that  fluid  element  travels  short  distance  within  short 
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time  interval,  and  that  the  information  carried  by  this  fluid  element  is  highly  preserved. 
The  displacement  vector  is  determined  by  finding  the  highest  correlation  between  two 
regions  on  two  successive  images  of  the  flow.  The  technique  works  as  follows.  First,  a 
small  window  called  the  source  window  is  selected  on  the  first  image.  Next,  this  source 
window  is  spatially  correlated  within  a  larger  window  on  the  second  image  called  the 
roam  window.  The  roam  window  is  centered  on  the  same  location  as  the  source  window 
as  shown  in  Figure  5. 


Roam  Window 
(Image  2,  t2=t1+Af) 

Soui  ce  Window 
(Image  1,  t=U) 


|  Displacement 
(Ax,  Ay) 


I _ I 


Figure  5.  The  source  and  the  roam  windows  of  the  image  correlation  velocimetry  technique. 

Next,  a  spatial  correlation  coefficient  C(Ax,Ay)  between  the  source  window  and  all  same 
size  windows  within  the  roam  window  is  calculated  as  a  function  of  the  pixel 
displacement  between  the  windows.  The  location  of  the  maximum  in  the  correlation 
coefficient  defines  the  displacement  vector  with  one  pixel  accuracy.  The  sub-pixel 
accuracy  could  be  achieved  via  two-dimensional  polynomial  fit  of  the  correlation 
coefficient  function  (Zheng  and  Klewicki,  2000).  The  equations  for  the  correlation 
coefficient  calculation  are  as  follows: 

C(A*,  Ay)  =  IsIr  ~  Is  Ir  (4) 

<7  <7 
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(5) 

(6) 


^MP-(i+ (7) 


Where  lr  and  Is  are  intensity  fields  of  source  window  and  the  displaced  window  inside  the 
roam  window,  and  M  x  N  is  the  size  of  the  source  window  (usually  square  window  is 
selected  so  that  M  -  N). 

Before  the  ICV  script  was  used  for  image  analysis  of  actual  microchannel  flows  it 
was  tested  with  artificially  designed  image  pairs.  The  first  image  in  the  test  pair  was  an 
image  acquired  with  our  micro  LEF  system  at  the  exit  of  the  microchannel.  The  second 
test  image  was  obtained  from  the  first  image  by  shifting  all  pixels  according  to  a 
prescribed  parabolic  velocity  profile.  The  ICV  script  was  tested  by  running  it  with 
different  sizes  of  the  source  window,  as  well  as,  analyzing  image  pairs  with  added 
random  noise. 


Figure  6  compares  the  profile  used  for  shifting  pixels  on  the  second  image  of  the 
test  pair  with  the  pixel  shift  profiles  obtained  by  ICV  script.  The  ICV  pixel  shift  profiles 
were  obtained  with  different  sizes  of  the  source  window.  It  is  seen  from  Figure  2  that  for 


Experimental  Study  and  Numerical  Modeling  of  Gas  Flow  in  Microchannels  and  Micronozzles 


20 


HCET-2003-M018-001-04 


all  used  source  window  sizes  pixel  shift  profiles  obtained  by  ICV  script  exactly  follow 
the  profile  used  for  creating  the  second  image. 

Figure  7  shows  velocity  vector  maps  for  the  test  image  pair  with  random  noise 
added  to  the  second  image.  Test  images  were  8-bit  gray  scale  images  with  256  gray 
levels.  Random  noise  of  5%,  7%,  and  10%  of  256  gray  scale  was  added  to  each  pixel  of 
the  second  image.  It  is  seen  from  Figure  3  that  5%  random  noise  does  not  significantly 
affect  the  vector  map  with  only  few  vectors  detected  incorrectly.  For  7%  noise  the  effect 
is  much  more  significant,  especially  for  small  pixel  shift  (low  velocity)  areas.  Finally  for 
10%  noise  the  vector  map  is  almost  completely  destroyed,  with  only  small  number  of 
vectors  detected  correctly. 

From  performed  ICV  script  testing  one  can  conclude  that  the  script  works  very 
accurately  and  even  could  tolerate  up  to  5%  random  noise  in  an  image  pair  without 
significant  deterioration  of  the  resulting  vector  map.  It  was  also  noted  that  the  effect  of 
adding  noise  was  more  pronounced  in  the  areas  of  small  pixel  shifts  (low  velocity). 
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Figure  6.  Comparison  of  the  shifting  profile  with  profiles  obtained  by  the  ICV  script. 
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Figure  7.  Velocity  vector  maps  for  test  image  pair  with  random  noise  added  to  the  second  image:  a) 
no  noise;  b)  5%  noise;  c)  7%  noise;  d)  10%  noise. 
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3.1.4  Flow  Velocity  Measurements  using  ICV  technique 

Upon  successful  testing  with  artificially  designed  image  pairs,  the  ICV  script  was 
used  for  velocity  field  measurements  for  free  jet  flow  exiting  from  the  microchannel  into 
air  under  normal  room  condition  (without  using  the  test  cell).  In  order  to  create  acetone 
concentration  gradient  for  velocity  detection  by  ICV  the  experimental  run  started  with 
nitrogen  gas  flowing  through  the  acetone  bubbler.  Than  the  flow  was  diverted  to  bypass 
the  bubbler  and  simultaneously  the  flow  rate  of  nitrogen  was  increased  to  supplement  for 
acetone  and  keep  total  flow  rate  constant.  Image  acquisition  was  initiated  simultaneously 
with  flow  diversion  and  a  number  of  image  pairs  were  acquired  capturing  the  process  of 
acetone  concentration  decreasing  from  its  original  value  to  zero.  Selected  pairs  were 
subsequently  analyzed  using  ICV  script.  Time  interval  between  the  images  in  one  pair 
was  selected  based  on  the  gas  flow  rate  and  desired  displacement.  For  example,  for 
nitrogen  flow  of  48  seem  average  expected  velocity  is  2.25  m/s  and  for  the  desired 
displacement  of  60  pixel  the  time  interval  between  the  images  in  the  pair  should  be  35  ps. 
Obtained  flow  fields  are  shown  in  Figure  8  to  Figure  1 1 .  Part  a)  of  each  figure  shows 
velocity  vectors  obtained  by  ICV  superimposed  on  the  LIF  signal  captured  on  the  first 
image  of  the  pair.  Shown  LIF  signals  are  pseudo  colored  proportional  to  the  strength  of 
the  LIF  signal  and  thus  acetone  concentration.  Part  b)  of  each  figure  shows  averaged 
vector  map.  On  these  figures  the  channel  exit  is  on  the  left  margin  of  the  figure  and  the 
gas  is  flowing  from  left  to  right.  Shown  field  of  view  is  0.825  mm  vertically  (slightly 
wider  than  channel  exit  width  of  0.5  mm)  and  1.7  mm  horizontally.  Time  interval 
between  the  images  in  one  pair  for  Figure  8  to  Figure  10  was  35  ps,  and  for  Figure  11  it 
was  8  ps.  To  obtain  flow  fields  on  Figure  8  to  Figure  11,  a  square  source  window  of  31 


Experimental  Study  and  Numerical  Modeling  of  Gas  Flow  in  Microchannels  and  Micronozzles 


23 


HCET-2003-M018-001-04 


pixel  side  was  used.  Since  flow  direction  is  known  from  the  experiment,  possible  source 
window  displacements  were  searched  in  the  flow  direction  within  rectangular  roam 
window  with  61  pixel  width.  Length  of  the  roam  window  was  selected  based  on 
maximum  expected  displacement  based  on  the  flow  conditions  and  time  between  image 
frames,  those  length  were  55  pixel,  90  pixel,  135  pixel,  and  70  pixel  for  Figure  8  to 
Figure  11  respectively.  Table  1  summarizes  average  measured  velocities  and  shows  their 
comparison  with  estimation  based  on  the  set  gas  mass  flow  rate. 


mm  pixel 


0 

0.12  50 
0.25  100 

0.37  150 

0.5  200 

0.62  250 
0.75  300 


100  2  0  0  300  4  0  0  500  600  pixel 

0  0.25  0.5  0.75  1  1.24  1.49  mm 


Figure  8.  Free  jet  flow  velocity  field  obtained  via  ICV  processing  of  acetone  LIF  images.  Expected 
average  velocity  0.75  m/s.  a)  Color  coded  acetone  LIF  image  is  superimposed  on  the  velocity  vector 
field,  b)  Averaged  vector  map. 
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Figure  9.  Free  jet  flow  velocity  field  obtained  via  ICV  processing  of  acetone  LIF  images.  Expected 
average  velocity  2.25  m/s.  a)  Color  coded  acetone  LIF  image  is  superimposed  on  the  velocity  vector 
field,  b)  Averaged  vector  map. 
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Figure  10.  Free  jet  flow  velocity  field  obtained  via  ICV  processing  of  acetone  L1F  images.  Expected 
average  velocity  3.6  m/s.  a)  Color  coded  acetone  LIF  image  is  superimposed  on  the  velocity  vector 
field,  b)  Averaged  vector  map. 
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Figure  11.  Free  jet  flow  velocity  field  obtained  via  ICV  processing  of  acetone  LIF  images.  Expected 
average  velocity  5  m/s.  a)  Color  coded  acetone  LIF  image  is  superimposed  on  the  velocity  vector  field, 
b)  Averaged  vector  map. 


Table  1.  Comparison  of  measured  flow  velocities  with  estimation  based  on  the  set  gas  mass  flow  rate. 


Gas  Mass  Flow 
Rate,  seem 

Estimated  average 
velocity,  m/s 

Measured  average 
velocity,  m/s 

Percent  deviation,  % 

16 

0.75 

0.66 

12.4 

48 

2.25 

2.05 

8.9 

77 

3.6 

3.16 

12.2 
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5 

4.42 

11.6 
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3.1.2  Molecular  Tagging  Velocimetry 

In  addition  to  velocity  measurements  by  ICV  technique,  molecular  tagging  velocimetry 
(MTV)  tests  were  conducted.  For  MTV  tests  following  arrangement  was  used.  The  laser 
beam  of  Nd:YAG  laser  was  narrowed  with  plano-convex  lens  from  4  mm  diameter  down 
to  about  100  pm  diameter.  The  laser  was  arranged  to  fire  close  to  the  microchannel  exit 
with  the  beam  perpendicular  to  the  flow  direction  as  shown  schematically  in  Figure  12. 
Nitrogen  gas  flow  was  constantly  seeded  with  acetone  vapor  to  produce  LEF  signal.  As 
illustrated  in  Figure  12,  sets  of  30  to  60  images  were  taken  simultaneously  with  laser 
shots  ( t-to )  and  with  a  time  delay  after  the  laser  shots  (t=to+At).  The  displacement  of  the 
line  written  by  the  laser  at  t=to  was  determined  by  comparing  the  LIF  signal  profiles  in 
the  flow  direction  on  averaged  images  taken  with  zero  delay  and  At  delay.  To  compare 
the  profiles  they  were  fitted  with  Gaussian  function  using  least-squares  fitting,  and  the 
displacement  was  obtained  from  the  position  of  the  peaks  of  the  obtained  fits.  Gaussian 
fitting  was  selected  because  it  describes  laser  energy  distribution  in  the  beam.  Gas 
velocity  was  obtained  from  the  measured  displacement  and  a  know  delay  for  the  images 
taken  after  the  laser  shot. 
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Figure  12.  Schematic  diagram  of  molecular  tagging  velocimetry  arrangement. 


Typical  line  written  into  the  flow  is  shown  in  Figure  13,  on  this  figure  the  channel  exit  is 
on  the  left  and  the  flow  is  from  left  to  right.  Image  on  Figure  13  is  taken  with  zero  delay 
after  the  laser  shot.  Typical  LIF  signal  profiles  of  written  line  are  shown  in  Figure  14. 
Two  profiles  are  shown  in  Figure  14,  one  with  zero  delay  and  the  other  with  80  ns  delay. 
Both  profiles  are  shown  with  corresponding  least-squares  fitting  with  Gaussian  function. 
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Figure  13.  Line  written  into  the  flow,  image  taken  with  zero  delay  after  the  laser  shot 
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Figure  14.  LIF  signal  profiles  with  Gaussian  fits  of  the  line  written  into  the  flow. 

A  one  dimensional  velocity  profile  obtained  by  MTV  technique  is  shown  in  Figure  15. 
This  profile  was  obtained  for  the  microchannel  inlet  pressure  of  92  torr  and  the  pressure 
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in  the  test  cell  of  62  torr  (ambient  pressure  at  the  exit).  Using  these  conditions  and 
channel  dimensions  in  the  solution  for  volume  flow  rate  for  incompressible  flow  in 
rectangular  duct  (White,  1974,  eq.  (10))  corresponding  average  flow  velocity  of  61  m/s 
was  obtained. 
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where  a  and  b  are  channel  height  and  width  (in  our  case  of  square  channel  a  =  b ). 

This  estimation  compares  very  favorably  with  average  velocity  of  60  m/s  from  the  data  in 

\ 

Figure  15. 


Figure  15.  Velocity  profile  near  microchannel  exit  (Pta  =  92  torr,  P3mb  =  62  torr). 


Fluorescence  intensity  of  initially  written  lines  as  a  function  of  pressure  in  the  test 
cell  is  shown  in  Figure  16.  It  is  seen  from  the  figure  that  LIF  intensity  is  approximately 
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linearly  dependent  on  pressure.  This  approximately  linear  dependence  was  previously 
observed  by  Lempert  et  al.  (2002)  for  pressure  range  up  to  9  torr.  Our  data  confirms  their 
result  and  shows  that  such  linear  dependence  is  observed  for  pressures  up  to  100  torr. 


Figure  16.  Fluorescence  intensity  of  initially  written  lines  as  a  function  of  pressure  in  the  test  cell. 


3.1.3  Gas  Mass  Flow  Rate  Through  Micronozzle 

Sandwich  type  design  micronozzles  with  rectangular  cross  section  were  obtained 
from  MEMS-exchange  manufacturing  company.  Photograph  of  one  of  the  manufactured 
nozzles  is  shown  in  Figure  17.  The  shown  micronozzle  has  throat  width  of  about  100  pm, 
which  is  somewhat  larger  than  the  design  width  of  80  pm. 

In  order  to  verify  the  integrity  of  the  micronozzle  and  gas  supply  line  mass  flow 
test  was  conducted.  A  simple  set-up  was  assembled  for  the  purpose  of  this  testing.  The 
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set-up  consisted  of  gas  cylinder  with  regulator,  pressure  transducer,  and  mass  flow  meter 
on  the  inlet  side  of  the  nozzle  and  the  other  mass  flow  meter  on  the  outlet  side  of  the 
nozzle.  Figure  18  shows  measured  flow  rates  before  and  after  the  micronozzle  as  a 
function  of  plenum  pressure.  It  is  seen  from  this  figure  that  for  pressures  up  to  40  PSI  gas 
flow  rate  before  and  after  the  micronozzle  are  the  same  within  the  measurement  error. 
Thus,  there  is  no  leaks  within  the  micronozzle  or  in  the  gas  supply  line. 


Figure  17.  Top  view  of  rectangular  cross  section  sandwich  type  micronozzle. 
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Figure  18.  Mass  flow  measurements  for  100  pm  throat  width  micronozzle. 


3.2  Numerical  Modeling 

3.2.1  Modeling  of  Gas  Flows  In  a  DeLaval  Nozzle  Using  Continuum  Model 
A  continuum  model  based  on  Navier-Stokes  equations  was  used  to  investigate  gas  flows 
in  a  DeLaval  nozzle.  Parametric  studies  included  the  effect  of  geometric  scaling  via 
scaling  down  of  throat  diameter  for  various  throat  Reynolds  numbers  and  the  effect  of 
chemical  propellants  on  the  integrated  performance  of  the  nozzle.  A  correlation  was 
developed  to  predict  the  specific  impulse  for  any  given  throat  diameter  and  throat 
Reynolds  number.  Detailed  description  of  the  model,  its  validation,  and  parametric  study 
results  are  presented  in  following  sections. 


Experimental  Study  and  Numerical  Modeling  of  Gas  Flow  in  Microchannels  and  Micronozzles 


34 


HCET-2003-M018-001-04 


3.2.1. 1  Numerical  Model  and  Its  Validation 

To  perform  numerical  modeling  of  gas  flows  in  a  DeLaval  nozzle  2-D  axisymmetric 
compressible  form  of  Navier-Stokes  equations  was  employed.  The  governing  euations 
are  presented  below. 


Continuity  equation: 
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Experimental  work  of  Jamison  et  al  (2003)  was  used  as  the  reference  problem  for  model 
verification  and  validation.  The  geometry  of  the  nozzle  used  in  simulations  is  the  same 


Experimental  Study  and  Numerical  Modeling  of  Gas  Flow  in  Microchannels  and  Micronozzles 


35 


HCET-2003-M018-001-04 


as  that  of  the  nozzle  used  by  Jamison  et  al  (2003).  As  shown  in  Figure  19  the  nozzle  has  a 
throat  diameter  of  1  mm  and  an  exit  diameter  of  7.9  mm,  which  gives  an  expansion  ratio 
of  e  =  62.41,  the  converging  and  diverging  sections  of  the  nozzle  have  cone  angles  of  30° 
and  20°,  respectively.  The  2D  axisymmetric  model  has  been  meshed  with  quadrilateral 
elements  as  shown  in  Figure  20. 

In  order  to  ensure  accurate  solutions  of  the  governing  equations,  a  grid  must  be 
sufficiently  fine  and  the  results  must  be  independent  of  grid.  A  grid  independence  test 
was  carried  out  before  the  actual  test  cases  were  run.  The  numbers  of  nodes  on  axial  and 
radial  directions  was  increased  to  generate  different  grids.  Grid  independency  is  defined 
as  a  stage  at  which,  the  results  remain  unchanged  with  the  increase  of  mesh  size.  To  find 
the  grid  independent  results,  several  cases  were  run  at  Re^oai  =  2.6  with  no-slip 
conditions.  The  criterion  for  the  grid  independency  in  the  present  study  was  the  variation 
of  the  mass  averaged  velocity  at  the  nozzle  exit  and  also  the  variation  of  velocity  at  a 
point  in  the  throat  region. 


Figure  19.  Nozzle  geometry. 
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Figure  20.  Nozzle  grid. 


A  240  x  400  grid  was  selected  as  the  finest  one  and  used  as  the  reference  to 
calculate  the  percentage  difference  of  velocity  obtained  from  the  grids  of  7  x  12  to  240  x 
400.  Figure  21  shows  the  grid  independence  study  results,  where  the  left  hand  scale 
denotes  the  velocity  (m/s),  while  the  scale  on  the  right  hand  side  represents  the 
percentage  of  the  velocity  deviation  with  the  fine  grid.  It  can  be  observed  from  Figure  21 
that  for  a  150  x  250  grid,  the  percentage  deviation  in  the  mass  averaged  exit  velocity  with 
fine  grid  was  about  1.4  %  and  for  the  percentage  deviation  for  the  velocity  at  throat,  y  = 
0.25  mm  was  0.1624  %.  The  150  x  250  grid  was  selected  as  the  optimum  grid  for  the 
simulations. 
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Figure  21.  Grid  independence  study  results. 

To  validate  our  numerical  model  we  compared  the  nozzle  exit  thrust  as  a  function  of 
Reynolds  number  with  experimental  measurements  of  Jamison  et  al  (2003).  The  nozzle 
exit  thrust  is  defined  as  the  sum  of  velocity  thrust  (MeVe)  and  pressure  thrust  (PeAe)  and 
is  given  by  the  equation  (16). 

Tn=MtVe  +  P'Ae  (16) 

Comparison  of  numerical  modeling  and  experimental  data  are  shown  in  Figure  22.  The 
experimental  thrust  calculations  consisted  of  velocity  thrust,  while  for  the  present 
numerical  study,  calculations  for  exit  thrust  were  attempted  with  and  without  the  pressure 
component  employing  both  slip  and  no-slip  boundary  conditions.  It  can  be  observed  from 
Figure  22  that  the  simulation  results  of  nozzle  exit  thrust,  which  included  pressure 
component,  are  significantly  larger  than  the  experimental  data  at  low  Reynolds  number, 
and  are  similar  at  increasing  Reynolds  number.  A  possible  reason  for  the  discrepancy  at 
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low  Reynolds  number  could  be  due  to  the  neglect  of  maximum  backpressures  and 
boundary  layer  growth. 


Figure  22.  Comparison  of  the  nozzle  exit  thrust  measured  experimentally  and  predicted  numerically. 

The  simulation  results  for  nozzle  exit  thrust  excluding  pressure  component  under 
predicted  the  experimental  data  for  both  slip  and  no-slip  flows  with  an  error  varying  from 
12  to  19  %  for  Rethroat  =  5.164  ~  22.692,  which  provides  reasonable  validation  of  the 
numerical  model.  It  was  observed  that  the  thrust  value  without  pressure  component  for 
no-slip  flow  was  slightly  higher  than  that  predicted  by  slip  flow  model  for  all  Rethroat  = 
0.164  ~  22.692.  However,  the  percentage  deviation  between  the  no-slip  flow  and  slip 
flow  was  very  small  between  0.02  to  4.66  %  for  Rethroat  =  0.164  ~  22.61.  Since  the 
experimental  data  provided  only  the  thrust  measurements,  which  gives  a  little  information 
about  the  physical  process  taking  place  inside  the  nozzle,  initial  simulations  were  carried 
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out  with  continuum  modeling  with  no-slip  boundary  conditions.  These  simulations  were 
used  to  predict  the  local  Knudsen  number,  which  indicated  that  the  flow  was  in  a 
transitional  regime  from  inlet  to  outlet  for  Rethroat  =  0.164  and  0.219.  For  Rethroat  =  0.3846 
to  22.692,  the  flow  is  in  a  slip  regime  for  nozzle  inlet,  while  for  the  throat,  the  flow 
continued  to  be  in  a  transitional  regime  until  Rethroat  =  2.637  and  a  continued  to  be  in  slip 
regime  for  Rethroat  =  22.692.  In  this  validation,  it  was  observed  that  the  model  predicted 
close  values  to  the  experimental  data  even  though  the  flow  was  in  transitional  regime  for 
certain  areas  inside  the  nozzle.  This  indicates  the  validity  of  the  model  even  for  the  flows 
in  transitional  flow  regime. 

3.2.1. 1  Geometric  Scaling 

Simulations  were  conducted  to  study  the  effect  of  the  geometric  scaling  of  the  nozzle  on 
the  integrated  performance  of  the  nozzle.  Figure  23  shows  the  scaling  down  of  the  nozzle 
schematically.  Parametric  studies  were  conducted  for  the  flows  ranging  for  Rethroat  equal 
to  5,  20,  50  and  100.  In  the  present  research,  all  the  dimensions  of  the  nozzle  were  scaled 
in  same  proportion  and  in  the  process,  the  throat  diameter  was  scaled  from  10,000 
microns  to  10  microns  while  maintaining  same  expansion  ratio,  converging  and  diverging 
angles  to  investigate  the  effects  of  nozzle  throat  diameter  with  the  increase  of  Rethroat- 
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Larger  D 


Figure  23.  Geometric  Scaling  of  the  Nozzle 

Figure  24  shows  the  variation  of  specific  impulse  Isp  as  a  function  of  the  throat  diameter. 
It  can  be  observed  from  the  Figure  24,  that  for  each  Rethroat,  with  the  decrease  of  Dt,  there 
was  one  fixed  Dt  (turning  point)  beyond  which  the  increase  in  the  value  of  Isp  is  relatively 
small.  There  is  a  shift  of  the  turning  point  with  Rethroat  and  was  marked  with  a  smooth 
curve  on  the  figure.  This  curve  divides  the  entire  plot  into  two  zones  known  as  scale 
sensitive  region  and  scale  numb  region.  The  variation  of  Isp  with  the  decrease  of  Dt  is 
significantly  large  in  the  scale  sensitive  region,  when  compared  to  scale  numb  region. 
One  possible  reason  for  the  turning  point  shift  could  be  due  to  the  effect  of  geometric 
dimension  on  Isp  at  low  Rethroat-  This  effect  subsided  with  the  increase  of  Rethroat- 
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Figure  24.  Variation  of  Isp  with  D,  for  Different  Re^o-,,  (Re  =  Reuu-oat). 


3.2.1.2  Correlation  for  the  specific  impulse 

A  correlation  was  developed  to  understand  the  variation  of  specific  impulse(s)  with  the 
variation  of  throat  diameter  and  throat  Reynolds  number  and  is  given  by  the  equation: 

lsp  =  C  Re 3  D,-2  04  (17) 

Where  C  can  be  approximated  to  29.54.  It  was  observed  that  Isp  is  varying  in  negative 
power  with  the  throat  diameter.  Figure  25  shows  the  correlation  plot  for  various  throat 
diameters  and  throat  Reynolds  numbers.  Using  this  correlation,  specific  impulse  of 
helium  propellant  can  be  approximated  for  any  given  throat  diameter  and  throat  Reynolds 
number  with  an  average  error  of  10  %. 
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Figure  25.  Correlation  for  the  specific  impulse  (Re  =  Rethroat) 


3.2.1.3  Effect  of  Chemical  Propellants 

Simulations  were  conducted  to  study  the  effect  of  the  chemical  propellants  on  the 
performance  of  the  nozzles  for  various  throat  diameters  of  10,  1  and  0.1  mm  respectively. 
These  throat  diameters  were  chosen  to  understand  the  effect  of  chemical  propulsion 
system  on  the  dimension  of  the  scale  (i.e.  micro  and  macro  scale).  The  propellants  were 
helium,  nitrogen,  argon  and  carbon  dioxide.  One  way  of  measuring  the  effectiveness  of  a 
chemical  propulsion  system  is  by  calculating  the  specific  impulse  of  the  ejected  gas. 
Mathematically,  it  can  be  expressed  as 

(18) 

8 

The  following  are  the  merits  for  adopting  specific  impulse  for  performance 
evaluation  of  the  nozzle.  First,  it  gives  us  a  quick  way  to  determine  the  thrust  of  a  rocket, 
if  the  mass  flow  rate  through  the  nozzle  is  known.  Second,  it  is  an  indication  of  engine 
efficiency.  Out  of  two  different  rocket  engines  that  have  different  values  of  specific 
impulse,  the  engine  with  the  higher  value  of  specific  impulse  is  more  efficient  because  it 
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produces  more  thrust  for  the  same  amount  of  propellant.  The  rocket  weight  will  define 
the  required  value  of  thrust.  Dividing  the  thrust  required  by  the  specific  impulse  will  tell 
us  how  much  weight  flow  of  propellants  the  engine  must  produce.  This  information 
determines  the  physical  size  of  the  engine. 

The  throat  Reynolds  number  was  varied  from  5  to  2000  for  four  different 
propellants.  Figure  26,  Figure  27,  and  Figure  28  show  the  variation  of  specific  impulse 
for  throat  diameters  0.1,1  and  10  mm  respectively.  It  can  be  inferred  that  helium  is  best 
propellant,  followed  by  nitrogen,  carbon  dioxide  and  argon.  It  was  observed  that  for 
throat  diameter  equal  to  0.1  mm,  there  was  not  much  difference  in  the  value  of  specific 
impulse  for  argon  and  carbon  dioxide  when  Repeat  is  below  250.  Above  Repeat  =  250, 
carbon  dioxide  dominated  argon.  Similarly  for  throat  diameter  equal  to  1  mm,  up  to 
Rethroat  =  250,  there  was  not  much  difference  between  argon  and  carbon  dioxide  and 
above  this  value  carbon  dioxide  dominated  argon.  Finally  for  Dt  =  10  mm,  argon 
produced  higher  specific  impulse  when  compared  to  carbon  dioxide  up  to  Rethroat  =  400, 
beyond  which,  carbon  dioxide  generated  higher  specific  impulse  than  argon.  It  was  also 
observed  that  there  was  no  significant  increase  in  the  value  of  Isp  when  Rethroat  is  greater 
than  750  for  all  propellants  for  the  three  throat  diameters  investigated. 


Experimental  Study  and  Numerical  Modeling  of  Gas  Flow  in  Microchannels  and  Micronozzles 


44 


HCET-2003-M018-001-04 


Experimental  Study  and  Numerical  Modeling  of  Gas  Flow  in  Microchannels  and  Micronozzles 


45 


HCET-2003-M018-001-04 


Figure  28.  Variation  of  1^  with  Reuiroat  for  D,  =  10  mm. 

Specific  impulse  calculated  at  Dt  =  10  mm  was  smaller  for  gases  like  nitrogen, 
argon  and  carbon  dioxide  when  compared  to  specific  impulse  calculated  at  Dt  =  1  mm. 
The  results  also  showed  that  there  was  no  significant  improvement  in  the  value  of 
specific  impulse  for  gases  like  nitrogen,  argon  and  carbon  dioxide  when  the  throat 
diameter  is  decreased  from  1mm  to  0.1  mm. 

Simulations  were  performed  at  a  mass  flow  rate  of  3.09  x  10'7  kg/s  to  study  the 
effect  of  gases  on  the  Isp  value.  The  results  are  shown  in  the  Table  2.  From  the  table  it  is 
very  clear  that  helium  dominated,  followed  by  nitrogen.  Argon  and  carbon  dioxide 
predicted  specific  impulse  almost  close  to  each  other  but  argon  predicted  a  little  higher 
value. 
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Table  2.  Specific  Impulse  at  a  given  mass  flow  rate  of  3.09  x  10'7  kg/s. 


Gas 

Isp(s) 

He 

72.61 

n2 

25.988 

At 

20.94 

co2 

20.86 

3.2.2  Modeling  of  Three-dimensional  Compressible  Gas  Flow  in  Microchannnels 
Using  Continuum  Model 

Nitrogen  gas  flow  in  long  microchannel  with  square  cross-section  was  simulated 
numerically  with  a  three-dimensional  continuum  model  with  slip  and  no-slip  boundary 
conditions.  The  governing  equations  of  the  model  were  solved  by  a  control  volume 
method.  Numerical  model  was  validated  with  available  experimental  and  numerical 
results.  For  incompressible  flow,  it  was  found  that  when  £>*  was  less  then  60  pm,  slip 
boundary  condition  must  be  applied.  Analytical  expression  for  normalized  friction 
coefficients,  C*lc ,  i.e.  ratio  of  fRe  (slip)  and  JRe  (no-slip),  was  developed  on  the  basis  of 
flow  behavior  for  incompressible  flow.  For  compressible  flow,  parametric  study  was 
conducted  for  Z>*  =  1pm,  UDh  =  200  and  with  varying  pressure  ratio  (PR  =  1.5  -  5.0).  It 
was  found  that  as  the  pressure  ratio  increased  from  1.5  to  5.0,  compressibility  effects 
increased  while  the  rarefaction  effects  started  diminishing.  The  slip  effects  also  played  an 
important  role  in  the  friction  characteristics  of  microchannel  flows.  Analytical  expression 
for  normalized  friction  coefficients,  C*c,  i.e.  ratio  of  fRe  (compressible)  and  fRe 
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(incompressible),  was  developed  on  the  basis  of  flow  behavior  for  compressible  flow. 
Comparative  study  of  two-dimensional  and  three-dimensional  flow  was  also  conducted 
and  it  was  shown  that  two-dimension  assumption  for  the  compressible  flow  was  not  valid 
since  it  gave  higher  flow  velocities,  15%  to  45%  higher  then  3-D  and  lower  friction 
factor,  7%  to  12%  lower  then  3-D. 

Nitrogen  gas  flow  in  3D  square  microchannel  was  studied  numerically  (Figure 
29).  For  incompressible  flow,  following  Chen’s  (2004)  validation  model  for  the  three- 
dimensional  steady  laminar  incompressible  flow  in  microchannel,  the  channel  outlet 
pressure  was  maintained  at  atmospheric  pressure  with  a  constant  pressure  ratio  of  1.0001. 
The  channel  length  over  hydraulic  diameter  ratio  was  kept  as  UDh  =  1000  &  200  with  the 
channel  hydraulic  diameter  varied  from  Dh  =  200  to  0.01  pm.  As  the  pressure  ratio  was 
very  small,  the  compressibility  effect  was  negligible.  For  the  compressibility  and 
rarefaction  effects  study  the  UDh  -  200,  Dh  =  1  pm  were  used.  The  channel  outlet 
pressure  was  maintained  at  atmospheric  pressure,  with  a  varying  pressure  ratio  from  1.5 
to  5.0,  in  these  high-pressure  ratios  the  compressibility  effects  were  crucial.  The  inlet  was 
located  at  x  =  0.  The  complete  physical  model  description  and  fluid  properties  are 
tabulated  in  Table  3. 


Figure  29.  MicroChannel  geometry. 


Experimental  Study  and  Numerical  Modeling  of  Gas  Flow  in  Microchannels  and  Micronozzles 


48 


HCET-2003-M018-001-04 


Table  3.  Physical  model  dimensions  and  fluid  properties. 


Parameter 

Range  or  Mean 

Working  Fluid 

Nitrogen 

Length,  L/Dh 

200- 1000  pm 

Width,  W  =  H 

0.01  -  200  pm 

Pressure  Ratio,  PR  =Pj„/Pout 

1.0001,  1.5  -  5.0 

Outlet  Pressure,  Pout 

101.325  kPa 

Inlet  Temperature,  Tj 

297  K 

Wall  Temperature,  Tw 

297  K 

Knudsen  Number,  Kn 

0.001  -  0.065 

Absolute  Viscosity,  p. 

1.663  x  1 0'5  Ns/m2 

Specific  Gas  Constant,  R 

296.7  J/kg  K 

Ratio  of  Specific  heat,  y 

1.4 

Molecular  mass 

28.013  kg/kmol 

The  simulated  3-D  gas  flow  and  heat  transfer  is  assumed  to  be  steady  and  laminar.  The 
continuity,  momentum  and  energy  equations  for  the  compressible  gas  flow  were 


expressed  as  follows: 
Continuity 


dpu  dpv  dpw  _ 
dx  dy  dz 


(19) 


Momentum 
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The  gas  velocity  and  temperature  at  slip  boundary  condition  were  governed  by 
equations  as  discussed  above.  All  walls  were  assumed  to  be  adiabatic  at  slip  and  no-slip 
conditions.  For  the  case  where  flow  was  considered  incompressible  the  density  term  was 
considered  constant.  Pressure  boundary  conditions  were  used  at  inlet  and  outlet 
boundaries  with  outlet  boundary  having  atmospheric  pressure  while  static  pressure  was 
specified  at  the  inlet. 

The  governing  Navier-Stokes  equations  were  solved  with  finite  volume  method. 
The  continuity  and  momentum  equations  were  solved  with  second  order  upwind  implicit 
scheme,  and  the  pressure-velocity  coupling  with  SIMPLEC  algorithm.  Energy  equation 
was  also  solved  by  second  order  upwind  implicit  scheme.  As  the  hydraulic  diameters  of 
the  simulated  channels  were  mostly  of  the  order  of  10'6  m,  the  scaled  convergence 
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criterion  was  set  four  orders  less  then  the  hydraulic  diameter  for  all  parameter 
calculations.  The  3-D  microchannel  was  divided  into  small  hexahedral  mesh  volumes  for 
numerical  calculations. 

Grid  independence  study  was  performed  to  ensure  properly  sufficient  number  of 
grids  had  been  generated  near  the  walls  as  shown  in  Table  4.  The  grid  of  size 
30x30x600  was  chosen  as  optimum  grid,  with  further  increase  of  grid  size  resulting  less 
than  2%  increase  in  average  outlet  velocity. 


Table  4.  Grid  independence  study. 


Grid  Size  V0u, 

Difference  Grid  Size 
(%) 

v0M 

Difference 

(%) 

10X10  X  600  0.689014 

5.90 

30X30X150 

0.7345011 

4.5 

20X20X600  0.7279103 

0.80 

30X30X300 

0.7342024 

2.3 

30  X  30  X  600  0.7335787 

0.60 

30X30X600 

0.7335787 

0.1 

40  X  40  X  600  0.7355692 

0.40 

30X30X900 

0.73374 

0.01 

50  X  50  X  600  0.7364773 

0.00 

30X30X1200 

0.73375 

0 

The  numerical  technique  was  validated  first  with  analytical,  experimental  and 
numerical  results  from  the  literature,  then  the  simulations  had  been  carried  out  for 
different  hydraulic  diameters  for  incompressible  gas  flow  to  study  the  effects  of  hydraulic 
diameter  on  gas  flow  properties  and  then  compressibility  and  rarefaction  effects  were 
studied  for  the  compressible  gas  flow,  and  finally  the  3-D  effects  were  studied. 
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3.2.2.1  Numerical  Technique  Validation 

To  validate  the  numerical  technique,  the  numerical  simulations  were  done  on 
conventional-size  channel  (Arkilic  et  al.,  1997).  The  channel  height  was  fixed  at  1  cm. 
The  channel  width  was  varied  to  obtain  different  aspect  ratios.  The  channel  length  was 
1000  times  the  hydraulic  diameter.  The  working  fluid  was  nitrogen.  The  channel  outlet 
pressure  was  atmospheric  pressure  and  the  pressure  ratio,  PR,  was  1.0001.  The  friction 
factor  for  the  flow  in  conventional  size  channels  was  calculated  for  different  aspect  ratios 
and  the  calculated  fRe  were  compared  with  the  analytical  solution  and  results  obtained  by 
Chen  (2004)  as  shown  below. 


Table  5.  Friction  factor  for  conventional  size  channel 


Re 

HAV  fRe 

fRe 

fRe 

(Numerical, 

(Analytical, 

(Numerical,  Present) 

Chen  (2004)) 

Chen  (2004)) 

116 

1  57.08 

56.91 

57.45 

212 

0.5  62.44 

62.19 

62.89 

Analytical  solution  for  fully  developed  incompressible  laminar  channel  flow  is 
given  below  (Chen  et  al.,  1998),  the  symbols  um,pm  are  the  average  velocity  and  density 

over  a  cross  section,  respectively. 

W2  dP 

u  = - 

m  Up  dx 

Analytical  solutions  for  mass  flow  rate  Q,  friction  factor/ and  Reynolds  number  Re 
were  calculated  using  the  analytical  um: 
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(25) 
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(28) 


The  flow  was  considered  incompressible  since  the  pressure  ratio  used  in  the 
numerical  simulations  was  very  small.  The  numerical  simulations  were  in  good 
agreement  with  the  numerical  simulations  reported  by  Chen  (2004)  and  with  the 
analytical  results  reported  by  Arkilic  et  al.  (1997).  The  agreement  was  within  the  1%. 

To  further  validate  the  numerical  model,  we  simulated  compressible  microchannel 
flows  and  compared  the  simulated  pressure  distribution  with  the  available  analytical 
solution  and  the  experimental  data  (Arkilic  et  al.,  1997).  The  dimensions  of  this 
microchannel  were  1.2x40x3000  pm.  To  comply  with  the  test  condition,  the  channel 
outlet  pressure  was  fixed  at  the  atmospheric  pressure.  The  pressure  ratio  was  2.708  with 
nitrogen  as  a  working  fluid.  Due  to  high  pressure  the  flow  was  considered  compressible 
and  the  slip  effects  were  also  included.  The  results  matched  well  with  the  experimental 
(Arkilic  et  al.,  1997)  and  numerical  (Chen,  2004)  results  thus  validating  our  numerical 
technique.  In  the  present  simulation  the  Reynolds  number  for  the  flow  in  Figure  30  was 
0.104. 
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Figure  30.  The  longitudinal  pressure  distribution  for  nitrogen  at  PR  =  2.708. 

3.2.2.2  Effect  of  Geometrical  Properties 

In  the  following  study  of  incompressible  gas  flow  in  the  square  microchannel,  the 
geometrical  properties  like  hydraulic  diameter  and  channel  length  were  varied  and  their 
effects  on  flow  properties  were  reported.  Hydraulic  diameter  is  equal  to  height  or  width 
of  square  microchannel.  The  Pressure  Ratio  (PR  =. Pin/Pout )  was  set  to  1.0001.  The  outlet 
pressure  was  fixed  as  atmospheric  pressure. 

In  our  study  by  Mo  et  al.  (2005),  the  square  microchannel  had  UDh  =  1000,  PR  = 
1.0001.  The  channel  hydraulic  diameter  was  decreased  from  Z>*  =  100  to  0.01  pm.  We 
reported  simulation  results  of  JRe  =  54.8,  which  differed  from  the  analytical  solution  by 
3.85%.  We  reported  that  the  slight  difference  could  be  resulted  from  the  neglect  of 
sidewall  effects  in  the  2-D  analytical  solutions. 
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Present  results  were  obtained  for  square  microchannel  with  L/Dh  =  200,  PR  = 
1.0001.  The  channel  hydraulic  diameter  in  this  case  was  also  decreased  from  A,  =  100  to 
0.01  Jim  to  investigate  the  effect  of  L/Dh- 

Mo  et  al.  (2005)  reported  that  when  the  gas  flow  was  in  slip  regime,  the  friction 
factor  would  be  reduced  due  to  the  slip  flow  effect  compared  with  no-slip  wall  condition, 
the  slip  model  calculated  fRe  started  to  decrease  at  Dh=  60  ftm  with  the  decrease  of  Dh, 
where  the  flow  was  still  in  continuum  flow  regime. 

Normalized  friction  coefficient,  C]c  =  /7?esiip//R<?NosiiP  was  calculated  for  the 
incompressible  gas  flow  here  from  Mo  et  al.  (2005).  As  shown  in  Figure  31  our  results 
matched  well  with  these  results.  We  can  conclude  that  for  the  incompressible  gas  flow 
reducing  the  L/Dh  from  1000  to  200  had  negligible  effect  on  normalized  friction 
coefficient.  So  to  save  expensive  computational  time  we  could  use  L/Dh  -  200  for  the 
developed  flow  calculation  in  microchannel.  We  could  also  conclude  that  as  the  hydraulic 
diameter  increased  the  effect  of  slip  boundary  conditions  reduced. 

An  analytical  expression  for  normalized  friction  coefficients  was  developed  on 
the  basis  of  behavior  for  incompressible  flow  as  shown  in  Figure  31. 


.  _  kDh 
ic  ~ 


Dh+c 


(29) 


Within  the  present  ranges,  the  values  of  constants  k  and  c  were,  k  ~  54.6  ±0.5  and  c 


0.6  ±0.05. 
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Figure  31.  C*lc  versus  hydraulic  diameter,  C]c  -fResuplfReno^v. 

3.2.2.3  Compressibility  and  Rarefaction  Effects 

To  study  the  compressibility  effects  for  slip  and  no-slip  boundary  conditions  on 
the  three-dimensional  microchannel  flows,  we  simulated  flow  in  square  microchannel 
with  flow  conditions  specified  in  Table  3.  The  channel  outlet  pressure  was  fixed  at 
atmospheric  pressure.  The  inlet  pressure  ratio  was  varied  from  1.5  to  5.0. 

The  pressure  drop  in  the  channel  takes  place  in  order  to  overcome  the  frictional 
forces  along  the  walls.  The  effect  of  slip  and  no-slip  boundary  conditions  on  the  pressure 
variation  (P/Pout  is  the  ratio  of  local  pressure  to  outlet  pressure)  along  the  centerline  of  the 
channel  is  shown  in  Figure  32  for  compressible  flow.  It  could  be  observed  that  with  the 
increase  in  pressure  ratio  the  divergence  in  the  pressure  distribution  between  slip  and  no¬ 
slip  wall  condition  increases.  The  maximum  percent  difference  in  the  pressure 
distribution  between  two  wall  condition  was  almost  6%  at  PR  =  5.0.  Linear  pressure 
distribution  for  incompressible  no-slip  flow  is  also  shown  for  comparison  in  Figure  32. 
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As  shown  in  Figure  32,  when  the  pressure  ratio  was  small,  PR  =  1.5,  the  compressibility 
effect  was  small,  and  pressure  distribution  was  almost  linear  for  all  the  flow  conditions. 
But  when  the  pressure  ratio  became  higher,  the  pressure  distribution  became  highly  non¬ 
linear  due  to  the  compressibility  effects.  So  we  could  conclude  that  the  nonlinearity  of  the 
pressure  distribution  and  corresponding  compressibility  effect  increased  with  the 
increasing  pressure  ratio  for  3-D  flow  conditions.  As  observed  in  Figure  32,  the  pressure 
distribution  for  the  no-slip  was  more  non-linear  when  compared  to  the  slip  flow.  The 
reason  could  be  that  the  slip  flow  reduced  the  wall  friction  and  therefore  diminishes  the 
non-linearity  of  the  pressure  distribution  compared  to  no-slip  flow. 


Figure  32.  The  longitudinal  pressure  distribution  at  various  pressure  ratios. 


As  the  pressure  ratio  (PR  =  Pin/P0ut)  increases  the  mass  flow  rate  increases 


nonlinearly  for  compressible  slip  and  no-slip  gas  flows  while  it  increases  linearly  for 


incompressible  no-slip  simulations  as  shown  in  Figure  33.  This  is  because  as  the  pressure 
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ratio  increases  the  compressibility  effect  becomes  more  dominant  and  thus  affects  the 
mass  flow  rate  more  substantially.  As  observed  in  the  Figure  33,  the  mass  flow  rate  was 
higher  with  the  slip  boundary  conditions  than  with  the  no-slip  boundary  conditions.  This 
was  the  result  of  the  lower  friction  resistance  in  the  microchannel. 


Figure  33.  Mass  flow  rates  at  various  pressure  ratios. 

Since  for  compressible  gas  flow  the  density  decreases  from  inlet  to  outlet  and  Knudsen 
number  varies  inversely  with  density.  Therefore  the  Knudsen  number  of  a  microchannel 
flow  increases  along  the  longitudinal  axis  and  reaches  the  maximum  at  the  outlet  as 
shown  in  Figure  34.  Figure  34  shows  the  local  Knudsen  number  variation  from  inlet  to 
outlet  for  two  pressure  ratios  for  compressible  flow.  The  behavior  was  similar  for  all 
other  pressure  ratios.  The  Knudsen  number  variation  from  inlet  to  outlet  for 
incompressible  gas  flow  was  very  minimal  and  thus  can  be  assumed  to  be  constant. 
Figure  35  shows  the  Knudsen  number  behavior  with  increasing  pressure  ratio.  The  mass 
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weighted  average  Knudsen  number  was  calculated  for  microchannel  at  each  pressure 
ratio.  It  could  be  seen  that  with  the  increasing  pressure  ratios  the  average  Knudsen 
number  was  reducing  for  the  compressible  case  while  for  the  incompressible  case  it  was 
almost  constant.  The  difference  between  the  average  Knudsen  number  for  slip  and  that 
for  no-slip  compressible  flow  was  negligible.  It  could  be  observed  that  as  the  pressure 
ratio  increased  from  1.5  to  5.0  the  average  Knudsen  number  was  almost  in  the  no-slip 
regime,  Kn  =  0.012. 


Figure  34.  Local  Knudsen  number  variation  along  the  channel  length. 
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Figure  35.  Average  Knudsen  number  with  varying  pressure  ratio. 

As  the  pressure  ratio  increased  the  average  mass  weighted  Mach  number  increased  as 
shown  in  Figure  36  compressible  no-slip  case  had  the  least  average  Mach  number  then 
compressible  slip  case.  At  relatively  low  PR  =  Pm/Pout,  average  Mach  number  for 
compressible  slip  flow  was  the  highest;  while  at  relatively  high  Pjn/P0ut,  the 
incompressible  no-slip  flow’s  average  Mach  number  was  the  highest.  The  reason  could 
be  that  as  the  pressure  ratio  increased  the  compressibility  effects  became  more  dominant. 
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Figure  36.  Average  Mach  number  with  varying  pressure  ratio. 

Normalized  friction  coefficient,  C'c=  fRecomp/fReincomp ,  was  calculated  for  varying 
pressure  ratio  as  shown  in  Figure  37.  The  incompressible  friction  factor  was  the  friction 
factor  for  incompressible  no-slip  boundary  conditions  fully  developed  flow,  which  was 
constant.  The  compressible  no-slip  case  had  higher  friction  coefficients  and  it  increased 
slightly  with  the  increasing  pressure  ratios.  The  compressible  slip  case  had  friction 
coefficients  lower  than  that  for  no-slip  flow.  These  results  showed  that  the  slip  effect 
indeed  reduced  the  wall  friction  significantly.  The  results  also  proved  that  the  fRe  value 
of  slip  flow  depend  strongly  on  the  pressure  ratio  (or  Reynolds  number)  especially  when 
Pin/Pout  was  low.  As  mentioned  earlier  the  slip  effects  decreased  with  increasing  pressure 
ratio  since  at  higher  pressure  ratio  the  gas  density  was  increased  and  the  mean  free  path 
of  gas  molecules  were  reduced.  The  Knudsen  number  and  the  slip  effect  diminish  with 
decreasing  mean  free  path. 
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Figure  37.  Cc  versus  varying  pressure  ratio,  Cc  =  fReComp/fReIncomp. 


An  analytical  expression  for  normalized  friction  coefficients  was  developed  on  the  basis 
of  behavior  for  compressible  slip  flow  as  shown  in  Figure  37. 


C*  =  kxe 


(30) 


Where  ki  and  k2  were  the  constants,  for  our  case  they  were  ki  -1.31  ±0.05  and  k2  -  - 
0.78  ±0.02. 


Figure  38  shows  the  streamwise  velocity  at  centerline  for  pressure  ratio  2.0-5.0,  at  the 
channel  exit,  for  slip  and  no-slip  flow  conditions  with  the  incompressible  no-slip 
condition  also  presented.  It  could  be  observed  in  Figure  38  (a)  and  (b)  that  compressible 
slip  flow  conditions  had  the  maximum  centerline  velocity  when  PR  =2.0  -  3.0.  But  as  the 
pressure  ratio  increased  from  3.0  to  5.0  the  incompressible  slip  flow  condition  had  the 
maximum  centerline  velocity  as  shown  in  Figure  38  (c)  and  (d).  As  discussed  earlier  this 
kind  of  behavior  could  be  due  to  the  diminishing  slip  effects  with  increasing  pressure 
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The  corresponding  velocity  solutions  at  the  channel  centerline  were  plotted  in 
Figure  39.  The  slip  condition  solutions  consistently  show  a  higher  magnitude  of  velocity 
than  no-slip  condition  due  to  lower  shear  stress  at  the  walls.  The  nonlinear  distribution  of 
the  velocity  was  dependent  on  the  pressure  distribution  and  hence  varies  for  different 
pressure  ratios.  As  the  pressure  ratios  increased,  the  non-linearity  of  the  centerline 
velocity  distribution  became  prominent. 


Figure  39.  Centerline  velocity  distribution  normalized  with  speed  of  sound  at  the  inlet  Cfa  for  various 
pressure  ratios. 

Figure  40  shows  the  flow  development  along  the  microchannel  for  PR  =  3.5.  As  could  be 
observed,  the  incompressible  flow  was  fully  developed  but  the  compressible  flow  was  not 
fully  developed  for  slip  and  no-slip  conditions.  Figure  40  shows  increase  in  streamwise 
velocity  with  the  corresponding  rise  in  wall  velocity  due  to  slip.  As  Kn  increased  due  to 
lower  density  downstream,  the  slip  affects  increased. 


Experimental  Study  and  Numerical  Modeling  of  Gas  Flow  in  Microchannels  and  Micronozzles 


64 


HCET-2003-M018-001-04 


Figure  40.  Slip  and  no-slip  U-velocity  at  different  cross-sections  along  the  y-direction. 


3.2.2.4  Three  Dimensional  (3-D)  Effects 

To  investigate  the  3-D  effects,  compressible  gas  flow  for  two-dimensional 
microchannel  was  simulated  to  compare  with  3-D  simulation.  The  geometric 
configurations  was  UDh  =  200,  Dh  =  1  pm. 

This  study  was  conducted  with  slip  flow  boundary  conditions.  Figure  41 
compares  the  velocity  solutions  along  the  channel  centerline.  As  can  be  observed  for 
exactly  the  same  conditions  the  U-velocity  component  always  had  higher  magnitude  for 
2-D  flow  then  3-D  flow.  The  difference  between  U-velocity  component  for  3-D  and  2-D 
flow  was  in  the  range  of  15%  to  20%. 

This  showed  for  compressible  flows  the  assumption  of  2-D  microchannel  was  not 
accurate  and  it  could  not  get  acceptable  results  as  expected.  Similar  behavior  was 
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observed  for  V-velocity  component  along  y-direction  at  the  channel  centerline  as  shown 
in  Figure  42.  The  difference  between  V-velocity  component  for  3-D  and  2-D  flow  was  in 
the  range  of  20%  to  25%. 

Figure  43  shows  the  W-velocity  component  at  channel  centerline  for  3-D  case  at 
the  exit  of  the  channel.  The  magnitude  of  W-velocity  increased  with  increasing  pressure 
ratio.  This  is  three-dimensional  effect,  which  was  neglected,  in  two-dimensional  study. 


Figure  41.  Centerline  U-velocity  distribution  for  compressible  gas  flow  with  slip  flow  boundary 
conditions. 
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Figure  44  compares  C*c  for  2D  and  3D  case.  As  could  be  seen  from  Figure  44,  the  two- 
dimensional  microchannel  always  gave  lower  C*c  and  thus  deviating  from  the  actual 

values.  This  validated  our  point  that  one  could  not  use  two-dimensional  assumption  to 
predict  the  compressible  flow  accurately.  The  difference  in  the  3-D  and  2-D  results  was 


in  the  range  of  7%  to  12%. 


—  —  Compressible,  Slip,  3-D 

1.0 

. * .  Incompressible,  No-Slip,  3-D 

- ♦ - Compressible,  Slip,  2-D 

1.2 

- 

1.1 

- 

4 

O 

°0.9 

"  — 

-  H 

0.8 

:  ^  - - 

- 

0.7 

r  ^ 

0.6 

- 

0.5 

-L  i.  1  .  ,  1_  .  I  ,  ,  ,  .  1  ,  .  .  .  _  1 

1  2  4  5 

h  out 

Figure  44.  Cc  versus  varying  pressure  ratio,  C'c  -  fReComp/fRelncomp. 

3.2.2.S  Effects  of  Aspect  Ratio 

To  study  the  effect  of  varying  channel  aspect  ratio  we  varied  the  height  and  width  of  the 
microchannels.  To  maintain  dynamic  similarity,  we  kept  the  hydraulic  diameter  fixed  as 
1/xm  as  shown  in  Table  4.  The  maximum  Knudsen  numbers  at  the  channel  outlet  were 
less  than  0.06.  In  the  following  simulation,  the  slip  effects  were  included.  The  channel 
outlet  pressure  was  fixed  at  atmospheric  pressure,  and  the  inlet  pressure  varied  from  2.0- 
5.0. 
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Table  6.  Microchannels  dimensions  used  for  aspect  ratio  effect  study. 


AR 

W(fim) 

H(iim) 

Dh 

1 

1 

1 

1 

0.75 

1.167 

0.875 

1 

0.5 

1.5 

0.75 

1 

0.25 

2.5 

0.625 

1 

0.20 

3 

0.60 

1 

Figure  45  depicts  the  calculated  normalized  friction  constant,  C*,  as  a  function  of 
pressure  ratio  for  different  aspect  ratio.  As  shown  in  Figure  45,  C*  increased  with  the 
decrease  in  aspect  ratio.  As  the  aspect  ratio  decreases  the  channel  height  becomes 
narrower,  which  gives  rise  to  a  wall  friction.  As  the  aspect  ratio  was  decreased  to  0.20, 
the  strongest  slip  effects  increased  the  C*  to  almost  twice  then  at  aspect  ratio  of  AR  =  1. 


Figure  45.  Normalized  Friction  constant  as  a  function  of  pressure  ratio  for  varying  aspect  ratio. 
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Figure  46  depicts  the  increase  in  mass  flow  rate  with  reducing  aspect  ratio.  As  shown  in 
Figure  47,  as  the  aspect  ratio  increased  the  mean  velocity  reduced,  this  is  due  to  increase 
in  the  wall  friction  effects. 


Figure  46.  Mass  flow  rate  as  function  of  pressure  ratio. 


Figure  47.  Mean  velocity  as  function  of  pressure  ratio. 
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The  investigation  of  effects  of  aspect  ratio  on  the  laminar  friction  constant  for 
compressible  gas  flows  would  be  incomplete  if  we  don’t  investigate  the  effects  on  flow 
properties  with  increase  in  aspect  ratio.  The  aspect  ratio  was  increased  from  1  to  5, 
keeping  the  hydraulic  diameter  constant  ( Dh  =  1pm).  As  shown  in  Figure  48,  the 
normalized  laminar  friction  constant  is  minimum  for  square  microchannel  and  it 
increases  as  the  aspect  ratio  decreases  or  increases  from  1.  Figure  49  shows  the  effect  of 
aspect  ratio  on  the  mass  flow  rate,  as  the  aspect  ratio  deviates  from  1  the  mass  flow  rate 
increases.  As  shown  in  figures  the  behavior  is  symmetrical. 


Figure  48.  The  C*  with  varying  aspect  ratio  at  PR  =  2.0 
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Figure  49.  Mass  flow  rate  with  varying  aspect  ratio. 

The  effect  of  aspect  ratio  was  further  investigated  at  the  various  pressure  ratios. 
As  can  be  seen  in  Figure  50,  as  the  pressure  ratio  increases  the  normalized  laminar 
friction  constant,  C*  increases,  this  is  due  to  increasing  compressibility  and  diminishing 
rarefaction  effects.  It  is  interesting  to  observe  that  the  behavior  of  the  friction  constant 
versus  aspect  ratio  is  symmetrical  and  the  trend  continues  for  varying  pressure  ratios.  A 
narrower  or  wider  channel  has  larger  friction  at  the  walls  and  therefore  higher  friction 
factor.  The  friction  characteristic  for  compressible  flow  in  microchannels  is  much  more 
complicated  than  that  for  incompressible  flow  in  conventional-size  channels.  For 
incompressible  conventional  channel  flow,  the  normalized  laminar  friction  constant  is 
only  a  function  of  aspect  ratio.  But  as  shown  in  Figure  50,  for  compressible  flow,  the 
normalized  laminar  friction  constant  is  a  function  of  Knudsen  number,  compressibility 
(pressure  ratio)  and  aspect  ratio  in  microchannels. 
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Figure  51  shows  the  mean  flow  velocity  as  a  function  of  aspect  ratio.  It  can  be 
observed  that  as  the  pressure  ratio  increases  the  mean  flow  velocity  increases  due  to 
compressibility  effects.  And,  at  aspect  ratio  of  1 .0,  the  channel  has  the  maximum  mean 
flow  velocity. 


Figure  50.  The  C*  with  varying  aspect  ratio  at  various  pressure  ratios. 


Figure  51.  Mass  flow  rate  with  varying  aspect  ratio  at  various  pressure  ratios. 
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3.2.3  DSMC  Simulations  of  Axisymmetric  Gas  Flow  in  Conical  Micronozzles 
The  direct  simulation  Monte  Carlo  (DSMC)  method  of  Bird  (1976)  is  a  statistical 
simulation  technique  that  models  a  real  flow  governed  by  Boltzmann  equation  with  a  set 
of  representative  molecules.  Current  applications  of  DSMC  are  limited  to  near  continuum 
and  rarefied  flows  due  to  its  intensive  computational  cost  (Chung  et  al.  1995).  2D  and 
3D  gas  flows  in  micronozzles  were  investigated  by  DSMC  to  aid  the  development  of 
microelectromechanical  system  (MEMS)  applications  in  micro-scale  thrusters,  etc. 
(Ivanov  et  al.  1999,  Markelov  and  Ivanov  2001,  Alexeenko  et  al.  2002,  Boyd  et  al.  1992 
and  1994),  where  Rothe’s  (1971)  experiment  was  frequently  used  as  the  reference 
problem.  Chung  et  al.  (1995),  studied  low-density  nozzle  flow  by  DSMC  and  continuum 
methods  based  on  Rothe’s  (1971)  experimental  data.  It  was  reported  that  the  continuum 
method  could  provide  good  results  inside  the  nozzle  in  terms  of  density  while  DSMC 
could  provide  good  results  in  both  density  and  rotational  temperature. 

The  present  study  is  focused  on  the  performance  evaluation  of  DSMC  and 
continuum  methods  for  the  simulations  of  normal-density  gas  flows  at  low  Reynolds 
number  in  a  conical  micronozzle  with  throat  diameter  of  1  mm  (in  contrast  to  Rothe’s 
(1971)  5  mm),  in  order  to  explore  the  advantages  of  using  the  two  approaches  in  different 
flow  regimes. 

Experimental  work  of  Jamison  et  al.  (2003)  is  used  as  the  reference  problem 
(Figure  52).  The  throat  diameter  is  Dt  =  1  mm  and  exit  diameter,  De  =  7.9  mm,  which 
gives  an  expansion  ratio  of  e=  62.41.  The  half  angles  are  30°  and  20°  in  the  converging 
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and  diverging  sections  respectively.  The  total  length  of  the  nozzle  is  L  =  10.7  mm. 
Helium  gas  is  used  as  the  working  fluid. 


Figure  52.  Nozzle  geometry. 

The  Reynolds  number  is  calculated  at  the  nozzle  throat  as: 


(31) 


Knudsen  number  is  calculated  at  the  nozzle  inlet: 


Kn=- 


D, 


(32) 


Where  X  is  molecular  mean  free  path  and  calculated  as  (Mohamed  1999): 


X  = 


kT 


-j2na2P 


(33) 


The  flow  regimes  are  categorized  by  Knudsen  number  as  (Mohamed  1999):  Kn  < 
10'3,  continuum  flow;  10'3<  Kn  <  10'\  slip  flow;  10'1  <  Kn  <  10,  transitional  flow;  Kn  > 
10,  free  molecular  flow.  Important  parameters  used  in  the  simulations  are  listed  in  Table 
7. 
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Table  7.  Computation  Parameters. 


Throat 

Nozzle  Inlet 

DSMC 

Re 

P(Pa) 

Kn 

Cell  Size,  m 

Timestep,  s 

0.16 

50.127 

0.1398 

1.12E-04 

8.28E-08 

0.22 

58.501 

0.1198 

9.62E-05 

7.10E-08 

0.38 

78.859 

0.0889 

7.14E-05 

5.26E-08 

0.55 

87.616 

0.0800 

6.42E-05 

4.74E-08 

0.77 

99.651 

0.0703 

5.65E-05 

4.17E-08 

1.04 

119.781 

0.0585 

4.70E-05 

3.47E-08 

2.64 

206.087 

0.0340 

2.73E-05 

2.01E-08 

5.16 

275.927 

0.0254 

2.04E-05 

1.50E-08 

7.58 

390.045 

0.0180 

1.44E-05 

1.06E-08 

10.93 

423.106 

0.0166 

1.33E-05 

9.81E-09 

14.45 

463.021 

0.0151 

1.22E-05 

8.97E-09 

18.08 

539.380 

0.0130 

1.04E-05 

7.70E-09 

22.69 

695.976 

0.0101 

8.09E-06 

5.97E-09 

Bird’s  DS2V  DSMC  codes  are  used  for  the  calculations.  Nozzle  geometry  is  divided  into 
a  set  of  cells  with  the  cell  size  calculated  as  one  third  of  the  molecular  mean  free  path. 
Timestep  is  calculated  as  one  third  of  the  molecular  mean  collision  time.  Initially,  the 
number  of  simulated  particles  in  the  computational  domain  is  assigned  arbitrarily.  The 
number  ratio  of  simulated  molecules  to  real  molecules  leads  to  the  statistical  scatter, 
following  the  Poisson  distribution,  with  the  standard  deviation  on  the  order  of  the  inverse 
square  root  of  the  sample  size.  The  simulation  fluctuations  may  become  unstable  in  some 
flow  configurations  due  to  random  walk  effect  in  one  or  more  variables.  One 
characteristic  of  a  random  walk  is  that  displacement  from  mean  position  or  value 
increases  with  square  root  of  time.  Random  walk  can  arise  whenever  one  of  the 
molecular  quantities  is  conserved  only  on  the  average,  rather  than  exactly,  in  any  of  the 
simulation  procedures. 

The  core  of  the  DSMC  algorithm  consists  of  four  primary  processes: 
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1.  Move  the  particles 

2.  Index  and  cross-reference  the  particles 

3.  Simulate  collisions 

4.  Sample  the  flow  field. 

These  procedures  are  uncoupled  during  each  time-step.  Of  primary  importance  is 
the  selection  of  a  timestep  that  is  less  than  the  mean  collision  time.  Simulation  results  are 
independent  of  the  timestep  increment  as  long  as  this  requirement  and  the  cell-size 
requirement  on  the  gradient  resolution  are  satisfied. 

The  DSMC  technique  is  explicit  and  time-marching,  so  that  it  always  produces  a 
flow  simulation  that  is  unsteady.  For  an  unsteady  flow  application,  an  ensemble  of  many 
independent  computations  may  be  assembled  and  averaged  to  obtain  final  results  with  an 
acceptable  statistical  accuracy.  An  ensemble  average  (the  instantaneous  average  over 
area  or  volume  elements  of  an  arbitrarily  large  group  of  similar  systems)  is  commonly 
used  to  present  unsteady  DSMC  results.  To  simulate  a  steady  problem,  each  independent 
computation  proceeds  until  a  steady  flow  is  established  at  a  sufficiently  large  time,  and 
the  desired  steady  result  is  a  time  average  of  all  values  calculated  after  reaching  the 
steady  state.  With  a  cumulative  sample,  the  statistical  fluctuations  declines  with  the 
square  root  of  the  sample  size  and  it  should  be  possible  to  achieve  any  desired  level  of 
accuracy  by  continuing  or  repeating  the  simulation  to  build  up  the  size  of  the  sample  to 
the  required  magnitude.  The  gradual  approach  of  a  sampled  quantity  to  its  correct  value 
is  sometimes  referred  to  as  ‘convergence’.  Figure  53  shows  the  judgment  of  convergence 
by  observing  the  smoothness  of  flow  stream  line. 
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b)  Possibly  converged. 


Figure  53.  Judgment  of  convergence  by  observations. 

As  shown  in  Table  7,  with  the  decrease  of  Knudsen  number,  DSMC  cell  size  and 
timestep  are  both  decreasing  due  to  the  decrease  of  the  molecular  mean  free  path.  In 
continuum  and  slip  flow  regimes,  the  fluid  is  much  denser  and  the  rarefaction  effect  is 
negligible,  where  DSMC  is  much  more  time  consuming  than  Navier-Stokes  (NS) 
method.  The  nozzle  exit  thrust  simulated  by  DSMC  and  NS  methods  show  good 
agreement  with  experimental  data  at  throat  Reynolds  number  larger  than  10  (Figure  54) 
where  the  flow  is  in  slip  regime.  With  the  decrease  of  throat  Reynolds  number,  the  fluid 
density  reduces  and  the  rarefaction  effect  becomes  important.  In  this  flow  regime,  NS 
simulation  is  heavily  influenced  by  the  back  flow  at  the  exit  and  the  exit  pressure  is 
limited  to  a  fixed  value  (1  Pa  in  this  report),  resulting  in  larger  predictions  of  exit  thrust. 
By  assuming  the  exit  pressure  is  negligible,  the  exit  thrust  calculated  by  NS  method 
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demonstrates  a  good  agreement  with  the  experimental  results.  While  the  DSMC  predicts 
larger  thrust  than  the  experimental  data  with  a  similar  variation  trend  with  the  throat 
Reynolds  number. 

As  shown  in  Figure  55,  at  throat  Reynolds  number  of  18.08,  the  pressure,  Mach 
number  and  temperature  contours  predicted  by  DSMC  and  NS  methods  are  similar  in 
pattern,  except  the  magnitude  and  smoothness  of  distributions.  The  exit  pressure  is  0  in 
DSMC  and  a  very  small  value  (close  to  1  Pa)  in  NS.  The  pressure  difference  is  not 
obvious  in  pressure  contour  (Figure  55  a),  but  it  causes  a  larger  Mach  number  predicted 
by  DSMC  (Figure  55  b).  The  temperature  distribution  near  the  wall  is  larger  in  DSMC 
predictions  (Figure  55  c). 
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Figure  54.  Nozzle  exit  thrust  vs.  throat  Reynolds  number. 
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a)  Pressure  contours  (Pa) 
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b)  Mach  number  contours 
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c)  Temperature  contours  (K) 


Figure  55.  Contours  predicted  by  DSMC  and  NS  methods  (throat  Reynolds  number  =  18.08). 


Figure  56  shows  the  effect  of  the  inlet  temperature  of  the  propellant  gas  at  Po  =  1 
Torr,  Dthroat=  1  mm.  With  higher  inlet  temperature,  the  molecular  kinetic  energy  is 
higher,  resulting  more  thrust  generated. 
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Figure  56.  Effect  of  the  inlet  temperature  of  the  propellant  gas  at  P0  =  1  Torr,  1  nun. 


3.2.4  Modeling  using  Unified  Flow  Solver 

This  section  contains  modeling  results  obtained  in  the  course  of  consulting  services  from 
CFD  Research  Corporation  (CFDRC)  and  presented  in  the  Final  Report  by  CFDRC 
entitled  ‘Testing  &  Demonstration  of  the  Unified  Flow  Solver  for  Micro  Flow 
Applications”.  The  team  has  been  working  with  CFDRC  in  selection  of  test  cases  and  in 
modeling  and  evaluation  of  the  Unified  Flow  Solver  (UFS)  for  gas  flows  in 
microchannels  and  micronozzles. 

The  Unified  Flow  Solver  is  developed  for  rarefied  and  continuum  flows.  The  UFS  code 
includes  deterministic  Boltzmann  solver  for  highly  non-equilibrium  flows  at  high 
Knudsen  number  and  continuum  (Euler  or  Navier-Stokes)  solvers  for  low  Knudsen 
numbers  (see  Figure  57).  Direct  Numerical  Solution  of  the  Boltzmann  equation  is 
preferable  to  DSMC  for  coupling  kinetic  and  continuum  models  and  for  simulations  of 
low  speed  flows. 
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Figure  57.  The  architecture  of  the  Unified  Flow  Solver. 

The  UFS  uses  continuum  breakdown  criteria  to  dynamically  select  kinetic  and  continuum 
domains  and  applies  appropriate  boundary  conditions  to  couple  the  Boltzmann  and 
continuum  solvers  at  the  boundaries. 

3.2.4.1  Boltzmann  and  Continuum  Solvers 

The  Boltzmann  transport  equation  (BTE)  for  a  particle  distribution  function /in  a 
one  component  atomic  gas  has  the  form: 

Yt  +  V(v/)  +  Vv-(a/)=  N-vf  (34) 

Here  r  is  spatial  coordinate,  v  is  velocity  vector,  and  a  is  acceleration  vector  The  left 
hand  side  of  Eq.  (34)  describes  the  free  motion  of  the  particles,  the  right  hand  side 
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describes  collisions  among  the  particles.  The  elastic  binary  collisions  are  described  by 
integral  operators 

y=j 

R  S 

J  J«r(U-5, |.®)|5-5, | mw.Warti.  (36) 

R  S 

Here  co  is  a  unit  vector  on  a  sphere  5  in  velocity  space,  and  dco  is  an  element  of  the  area 
of  the  surface  of  this  sphere,  where  (^',  £*')  represent  the  post  collisional  velocities 
associated  with  the  precollisional  velocities  ( £ ,  %  *)  and  the  collisional  parameter  a) ,  and 
the  kernel  <y  describes  details  of  the  binary  interactions. 

The  kinetic  gasdynamics  scheme  for  the  continuum  equations  is  based  on  finite 
volume  evaluation  of  semi-fluxes  at  cell  faces  using  analytical  values  for  a  Maxwellian 
distribution  function.  For  the  Navier-Stokes  equations,  the  semi-flux  scheme  is  based  on 
BGK  equation.  Directional  splitting  method  is  used,  and  flux  in  one  direction  is  evaluated 
by  means  of  the  construction  the  distribution  function  f0(t,x,£),  and  Maxvellian  g(x,t) 
around  the  cell  face  at  each  time  step.  For  a  ID  case  this  means: 

/0W  =  g‘[l  +  a‘x  -  T(a‘€  +  A1  )](1  -  H[x])  +  gr[  1  +  arx  -  r(ar<f  +  Ar)]H[x] , 

g(x)  =  So  [1  +  (1  -  H[x])dlx  +  H[x]drx  +  At] ,  (37) 

where  H[x]  is  the  Heaviside  step  function,  gl  and  grare  Maxwellian  distributions  at  the 

left  and  right  sides  of  the  cell  face,  g0  is  a  Maxwellian  at  the  cell  face,  and  functions 

al(r\al(r) ,  Al(r)  and  A  are  polynomials  with  local  coefficients.  The  fluxes  at  cell  faces  are 
computed  using  general  solution  of  the  BGK  equation  for  the  distribution  function /.  Note 
that  term  r(a,(r)£  +  Al(r))gUr)  defines  deviation  of  the  distribution  function  from 
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equilibrium,  where  T-filp.  This  procedure  gives  equations  for  mass,  momentum  and 
energy  transport  and  approximates  NS  equations.  If  the  parameter  r  is  set  to  zero,  the 
Euler  kinetic  scheme  is  obtained. 

The  automatic  switch  between  the  Boltzmann  solver  and  the  continuum  solvers  is 
performed  using  binary  switches  that  belong  directly  to  the  “cell”  structure  and  are 
accessible  everywhere  throughout  the  UFS  code.  The  memory  for  the  Boltzmann  solver 
is  allocated  only  in  the  kinetic  regions  of  the  computational  domain.  The  rest  of  the 
computational  domain  is  solved  by  the  continuum  (Euler  or  NS)  solvers.  The  switch  is 
based  on  gasdynamic  macroparameters  and  so  is  known  everywhere  in  the  computational 
domain. 

The  following  switching  criteria  are  available  in  the  source  code: 


where  pM'PyyPu  are  appropriate  components  of  the  non-equilibrium  stress  tensor,  p  is 
the  pressure,  p  is  density,  T  is  temperature,  u,  v,  w  are  appropriate  component  of  velocity 
(all  values  are  given  in  dimensionless  form).  If  S  is  greater  then  a  threshold  value,  then 
the  kinetic  solver  must  be  used.  The  applicability  of  different  criteria  and  the  ways  to 
choose  the  threshold  value  is  currently  being  studied.  So  far,  it  was  obtained 
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experimentally  that  criterion  SNS  gives  correct  kinetic  regions  and  allows  successfully 

couple  NS  kinetic  solver  with  Boltzmann  solver. 

In  order  to  match  the  Boltzmann  and  the  Euler  solutions,  a  special  technique  is 
used  which  consists  in  the  following.  In  a  narrow  layer  around  the  Boltzmann  region, 
Maxwellian  distribution  is  calculated  based  on  local  gasdynamics  macroparameters.  The 
advection  part  of  the  Boltzmann  solver  is  then  run  using  these  Maxwellian  distributions 
in  the  surrounding  layer. 

3.2.4.2  Simulation  of  MicroChannel  Flow 

Simulation  conditions  corresponding  to  experiments  of  Hsieh  et  al.  (2004)  were 
used  for  microchannel  flow  test  case. 

In  order  to  reduce  the  computational  effort,  the  simulation  domain  in  the  y-z  slice 
includes  only  a  quarter  of  the  domain  under  the  study  as  shown  in  Figure  58.  Two 
boundaries  in  the  y-z  slice  are  considered  as  symmetry  and  the  other  two  as  Neumann 
boundaries  on  which  diffusive  reflection  with  T  =  1  is  set  (see  Figure  58). 

The  inlet  and  exit  boundaries  are  treated  as  Inflow  and  Outflow  boundaries.  On 
the  Inflow  boundary,  the  density  and  pressure  is  set;  the  velocity  field  is  taken  from  the 
simulation  domain  values  so  that  normal  derivatives  of  the  velocities  are  zero.  On  the 
Outflow  boundary,  only  the  pressure  is  set;  the  velocity  and  density  field  is  taken  from 
the  simulation  domain  values  so  that  normal  derivatives  of  the  density  and  velocities  are 
zero.  This  is  a  usual  way  of  setting  boundary  conditions  for  simulations  of  flow  in  long 
tubes. 
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Results  of  the  simulation  are  shown  in  Figure  59  as  contour  plots  of  the  density, 
longitudinal  velocity,  and  temperature. 


Dlff.  reflection  T  =  1 


z 


Figure  58.  Simulation  domain  and  boundary  conditions  in  the  y-z  slice  of  the  studied  geometry. 
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Figure  59.  Results  of  microchannel  flow  simulation  by  Unified  Flow  Solver:  a)  density;  b) 
longitudinal  velocity;  c)  temperature. 


3.2.4.3  Simulation  of  Micronozzle  Flow 

Experimental  and  modeling  results  of  Bayt  (1999)  were  used  as  a  reference  case 
for  testing  UFS  for  micronozzle  flow  applications.  Following  boundary  conditions  were 
used  rho  =  1,  u  =  0.2*sqrt(5./6.*1.2)  and  T=  1.2  at  the  entrance  and  rho  =  0.05,  u  =  1.73 
and  T  =  0.2  at  the  exit.  Thus,  at  the  entrance,  the  Mach  number  is  set  to  0.2.  These 
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boundary  conditions  are  close  to  those  in  the  Bayt  thesis  (1999).  The  results  of 
simulations  using  the  gasdynamic  solver  are  shown  in  Figure  60.  For  comparison,  we 
show  the  Mach  number  profiles  from  the  Bayt  theses  in  Figure  61.  One  can  see  a  good 


comparison  with  the  results  of  the  UFS  code. 


c) 

Figure  60.  Results  of  simulation  of  the  Bayt’s  nozzle  using  the  Euler  solver:  a)  Density  contours;  b) 


Mach  number  contours;  and  c)  Temperature  contours. 
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Figure  61.  Mach  number  contours  for  the  16.9:1  nozzle  at  a  Reynolds  number  of  1914  (from  Bayt 
(1999)). 


The  above  simulation  was  repeated  with  Kn=0.001  to  demonstrate  the  automatic 
switch  between  the  gasdynamic  and  Boltzmann  solvers.  The  results  are  presented  in 
Figure  62.  One  can  see  that  there  is  a  considerably  large  kinetic  region  (shown  in  red  in 
Figure  62  a))  in  the  portion  of  the  nozzle  with  low  density. 
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Figure  62.  Results  of  micronozzle  flow  simulation  by  Unified  Flow  Solver  (Kn=0.001):  a)  kinetic  flag; 


b)  density  contours;  c)  Mach  number  contours;  d)  temperature  contours. 
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3.3  Summary  and  Conclusions 

A  test  set-up  capable  of  establishing  controlled  gas  flow  in  microchannels  and 
micronozzles  has  been  designed,  assembled,  and  tested.  Validation  of  the  set-up  was 
performed  by  conducting  flow  tests  with  pure  nitrogen  gas  and  comparing  measurements 
with  classical  correlations  for  fluid  flows  in  square  duct.  Reasonable  agreement  of 
experimental  data  with  correlations  was  observed  when  hydraulic  diameter  of  500  pm  to 
550  pm  was  used  in  the  correlations.  A  micro  gas  Laser  Induced  Fluorescence  system 
was  setup.  Image  analysis  script  utilizing  Image  Correlation  Velocimetry  (ICV) 
principles  was  written  and  tested  with  specially  designed  test  images  as  well  as  with 
images  with  added  random  noise.  The  velocity  profile  obtained  by  ICV  script  was  found 
to  follow  closely  the  test  profile  used  to  obtain  test  images.  The  script  was  capable  to 
tolerate  random  noise  of  up  to  5%  of  total  gray  scale.  Gas  velocity  measurements  in 
microscale  were  conducted  using  ICV  technique  for  free  jet  flow  exiting  from  the 
microchannel.  Maximum  deviation  of  experimentally  measured  average  velocity  via  LIF- 
ICV  technique  from  estimation  based  on  the  set  gas  mass  flow  rate  was  12.4%.  Velocity 
measurements  were  also  preformed  via  Molecular  Tagging  Velocimetry  technique  and 
produced  good  agreement  with  estimation  of  average  velocity  for  incompressible  flow  in 
a  rectangular  duct. 

Continuum  modeling  on  nozzle  flows  was  conducted  and  validated  by  available 
experimental  data.  A  numerical  study  has  been  performed  on  the  effects  of  throat 
diameter  via  geometric  scaling  down  of  throat  diameter  from  10  mm  to  0.1  mm  for  four 
different  throat  Reynolds  numbers  with  helium  as  test  propellant.  Variation  of  Isp  was 
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found  to  be  negligible  from  1mm  to  0.1  mm  when  compared  to  10  mm  to  1mm  of  throat 

diameter  for  a  given  throat  Reynolds  number  for  helium  flows.  Similar  behavior  for  Isp 

was  observed  for  other  propellants.  A  correlation  for  specific  impulse  involving  throat 

diameter  and  throat  Reynolds  number  was  developed  and  presented.  The  correlation  is 

capable  of  predicting  Isp  with  an  error  of  up  to  10  %  of  the  actual  value. 

In  a  numerical  study  of  gas  propellants  is  was  found  that  helium  proved  to  be  best 

among  nitrogen,  argon  and  carbon  dioxide  for  throat  Reynolds  number  varying  from  5  to 

2000  for  the  three  throat  diameters  investigated.  Nitrogen  followed  helium.  The 

comparative  performance  of  argon  and  carbon  dioxide  varied  depending  on  the  throat 

diameter  and  the  throat  Reynolds  number. 

In  the  modeling  of  three-dimensional  gas  flow  in  microchannnels,  for 
incompressible  model,  the  effect  of  variation  of  UDh  on  fRe  or  C)c  was  found  negligible. 

As  the  hydraulic  diameter  was  varied  from,  £>/,  =  0.01  to  100  pm  friction  factor  increased. 
An  analytical  correlation  was  developed  to  predict  the  normalized  friction  factor  variation 
with  change  in  pressure  ratio  for  incompressible  flow. 

Compressibility  and  rarefaction  effects  in  three-dimensional  microchannel  gas 
flow  were  studied  with  Reynolds  number  varying  from  0.0001  to  1.2,  pressure  ratio 
varying  from  1.5  to  5.0,  and  Knudsen  number  varying  from  0.001  to  0.06.  The  findings 
are  summarized  as  follows: 

•  Due  to  compressibility  effects,  the  pressure  drop  exhibits  an  unusual  nonlinear 
behavior  as  compared  to  that  in  a  larger  channel.  Rarefactions  effects  were 
important  but  as  the  pressure  ratio  (or  Reynolds  number)  increased  the 
compressibility  effects  became  major  dominant  factor. 

•  Knudsen  number  increased  from  inlet  to  outlet  for  compressible  flow  but 
Knudsen  number  decreased  as  the  pressure  ratio  increases. 
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•  With  increasing  pressure  ratio  the  normalized  friction  factor  increased  for  both 
slip  and  no-slip  flows.  For  the  studied  range  of  pressure  ratios,  the  fRe  for 
compressible  slip  flow  is  less  than  fRe  for  incompressible  case. 

•  An  analytical  correlation  was  developed  to  predict  the  normalized  friction  factor 
variation  with  change  in  pressure  ratio  for  compressible  slip  flow. 

•  Two  dimensional  assumption  was  not  valid  for  compressible  flow  study  in 
microchannel  since  it  gave  much  higher  velocities,  15%  -45%  higher  then  3-D 
and  lower  normalized  friction  factor,  7%  -12%  lower  then  3-D. 

Effects  of  aspect  ratio  on  compressible  gas  flow  in  microchannels  were  studied 

for  the  ranges  of  pressure  ratio  from  2.0  to  5.0,  Reynolds  number  from  0.01  to  1.2,  and 
aspect  ratio  from  0.2  to  5.0.  It  was  observed  that  the  change  from  square  to  rectangular 
microchannel  results  in  the  increase  of  the  normalized  friction  constant.  The  increase 
could  be  almost  two  times  higher  for  narrower  rectangular  microchannels  compared  to 
the  square  microchannel.  For  compressible  flow,  the  normalized  laminar  friction 
constant  was  observed  to  depend  on  Knudsen  number,  compressibility  (pressure  ratio) 
and  aspect  ratio. 

In  DSMC  simulations  of  axi  symmetric  gas  flow  in  conical  micronozzles  it  was 
found  that  at  throat  Reynolds  numbers  larger  than  10,  DSMC  and  continuum  methods 
predicted  closer  nozzle  exit  thrusts  to  the  experimental  data.  Pressure,  Mach  number  and 
temperature  contours  were  predicted  by  DSMC  and  continuum  methods  to  be  similar  in 
pattern,  but  different  the  magnitude  and  smoothness  of  distributions.  DSMC  predicted 
larger  Mach  number  and  temperature  near  the  wall  and  closer  results  at  throat  Reynolds 
numbers  less  than  10  where  the  flow  is  switching  from  slip  to  transitional  regimes.  A 
study  of  propellant  gas  temperature  effect  showed  that  at  Po  =  1  Torr,  Dthroat=  I  mm,  the 
higher  the  gas  temperature,  the  larger  the  thrust  is  generated. 
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A  Unified  Flow  Solver  that  utilizes  hybrid  approach  using  deterministic 
Boltzmann  solver  for  highly  non-equilibrium  flows  at  high  Knudsen  number  and 
continuum  (Euler  or  Navier-Stokes)  solvers  for  low  Knudsen  numbers  was  tested  and 
demonstrated  for  gas  flows  in  microscale.  Initial  examination  showed  that  a  reasonable 
agreement  with  the  data  from  open  literature  was  obtained  for  the  micronozzle  test  case. 
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4.0  Nomenclature 


i 


A*  Exit  Area  (m2) 

D,  Throat  Diameter  (mm) 

Dh  Hydraulic  Diameter  (mm) 

De  Exit  Diameter  (mm) 

g  Acceleration  of  gravity  (9.8  m/s2) 

Isp  Specific  Impulse  (s'1) 

k  Boltzmann’s  Constant  (1.38*10‘23  J/K) 

K  Thermal  Conductivity  (W/m-K) 

Kn  Knudsen  Number 

L  Characteristic  Length  (m) 

Me  Mass  Flow  Rate  (Kg/s) 

P  Pressure  (Pa) 

Q  Volume  flow  rate  (m3/s) 

Rethroa,  Throat  Reynolds  Number 
T  Temperature  (K) 

Tn  Nozzle  Exit  Thrust  (N) 

p  Fluid  Density(Kg/m3) 

X  Mean  Free  Path  (m) 

p  Molecular  Viscosity  (Kg/m-s) 

ctc  Collision  Diameter  of  the  Molecules  (m) 
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6.0  Appendix 


Table  8.  Experimental  data  and  correlation  predictions  for  N2  flow  in  a  square  microchannel  (see 
Figure  4). 


Pin,  torr 

Pom,  torr 

Q,  seem 

dP,  torr 

dP corr  Dh=500«m,  tOIT 

dPcorr  Dh=550wm,  tOIT 

16.80556 

12.13773 

14.68102 

4.66782 

9.067 

6.193 

27.717 

20.87664 

26.96391 

6.84036 

10.38 

7.09 

36.80702 

28.36647 

39.50738 

8.44056 

11.771 

8.04 

44.55369 

34.73196 

50.86811 

9.82173 

12.827 

8.761 

51.96964 

40.94871 

62.70616 

11.02093 

13.893 

9.489 

59.20849 

46.93731 

76.39973 

12.27118 

15.275 

10.433 

66.29013 

52.94336 

88.76407 

13.34678 

16.243 

11.094 

73.11198 

58.62907 

100.6728 

14.48291 

17.091 

11.673 

79.49758 

64.0252 

113.57422 

15.47238 

18.134 

12.386 

85.92691 

69.4262 

125.21704 

16.50071 

18.497 

12.634 

92.1764 

74.68936 

137.02256 

17.48704 

18.868 

12.887 

98.85651 

80.28471 

148.65989 

18.5718 

19.088 

13.037 

Table  9.  Velocity  profile  near  microchannel  exit  obtained  by  MTV  technique  (Pin  =  92  torr,  Pamb  =  62 
torr,  see  Figure  15). 


Distance,  /tin 

V,  m/s 

-292 

80.92 

-219 

62.24 

-146 

37.35 

-73 

59.13 

0 

96.48 

73 

93.36 

146 

46.68 

219 

24.9 

292 

80.92 
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Table  10.  Fluorescence  intensity  of  initially  written  lines  as  a  function  of  pressure  in  the  test  cell  (See 
Figure  16). 


Pceii,  torr 

LIF  intensity 

12.1 

2 

23.7 

3 

35.7 

5 

50.2 

8 

67.7 

11 

75.9 

12 

96.2 

15 

Table  11.  N2  mass  flow  rate  through  the  100  pm  throat  width  micronozzle  (see  Figure  18). 


Pin,  psi 

Qm  in,  kg/s 

Qmouti  kg/s 

Error  Qm  in,  kg/s 

Error  Qra  out,  kg/s 

2.80666 

6.929E-8 

8.333E-8 

1.942E-8 

1.942E-8 

8.79558 

1.84E-7 

1.667E-7 

1.942E-8 

1.942E-8 

12.00444 

2.43E-7 

2.083E-7 

1.942E-8 

1.942E-8 

14.87842 

2.697E-7 

2.5E-7 

1.942E-8 

1.942E-8 

18.11976 

3.201E-7 

2.917E-7 

1.942E-8 

1.942E-8 

21.1042 

3.544E-7 

3.333E-7 

1.942E-8 

1.942E-8 

24.00094 

3.883E-7 

3.542E-7 

1.942E-8 

1.942E-8 

26.89764 

4.12E-7 

3.958E-7 

1.942E-8 

1.942E-8 

30.05596 

4.735E-7 

4.375E-7 

1.942E-8 

1.942E-8 

33.02218 

5.144E-7 

4.583E-7 

1.942E-8 

1.942E-8 

35.9765 

5.276E-7 

5E-7 

1.942E-8 

1.942E-8 

38.83834 

5.772E-7 

5.417E-7 

1.942E-8 

1.942E-8 
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Table  12.  Nozzle  exit  thrust:  comparison  of  experimental  data  with  continuum  Navier-Stokes 
solution  with  slip  and  no-slip  boundary  conditions  (see  Figure  22). 


Exit  Thrust  (Tn),  N 


Expt 


0.1648353 


0.2197804 


0.3846155 


0.5494507 


0.7692309 


1.043956 


2.637363 


5.164835 


7.582418 


10.93407 


14.45055 


18.07692 


22.69231 


^ [Till 


4.54E-07 


8.82E-07 


2.01  E-06 


2.76E-06 


4.45E-06 


6.45E-06 


1.72E-05 


3.92E-05 


6.67E-05 


1.13E-04 


1.69E-04 


2.32E-04 


3.19E-04 


1.06389E-07 


1.79329E-07 


4.85461E-07 


8.97384E-07 


1.57842E-06 


2.6 1363 E-06 


1.14733E-05 


3.24186E-05 


5.73122E-05 


9.69756E-05 


1.44584E-04 


1.94417E-04 


2.57116E-04 


1.12807E-07 


1.90828E-07 


5.10002E-07 


9.37582E-07 


1.64526E-06 


2.71020E-06 


1.18251E-05 


3.29897E-05 


5.82437E-05 


9.91597E-05 


1.48100E-04 


1.96718E-04 


2.60746E-04 


4.91254E-05 


4.91987E-05 


4.95075E-05 


4.99219E-05 


5.06056E-05 


5.16453E-05 


6.05172E-05 


8.1481  IE-05 


1.06387E-04 


1.46046E-04 


1.93658E-04 


2.44574E-04 


3.35296E-04 


4.91318E-05 


4.921 10E-05 


4.9533  IE-05 


4.99633E-05 


5.06739E-05 


5.17421E-05 


6.08740E-05 


8.20540E-05 


1.07319E-04 


1.48246E-04 


1.97190E-04 


2.68739E-04 


3.76755E-04 


Table  13.  The  longitudinal  pressure  distribution  for  nitrogen  at  PR  =  2.708  (See  Figure  30). 


Analyt.  Incompr. 

(Arkillic,  1997) 


x/L  P ZP, out 


6.52E -09  2.7142 


1  I  1 


x/L 


6.52E-09 


0.1344262 


0.3704918 


0.6 


0.8295082 


P/Pout 


2.7142 


2.584232 


2.262435 


1.860188 


1.433188 


Numerical  Compress.  (Chen, 
2004) 

x/L 

P/Pout 

6.52E-09 

2.7142 

0.09836067 

2.578043 

0.1344262 

2.534725 

0.1967213 

2.460464 

0.2327869 

2.404768 

0.2655738 

2.355261 

0.2983607 

2.31813 

0.3344262 

2.262435 

0.4 

2.169609 

0.4360656 

2.113913 

0.4655738 

2.064406 

0.5016394 

2.00871 

0.5311475 

1.959203 

0.6327869 

1.785928 

0.6655738 

1.724044 

0.6983607 

1.668348 

0.7344263 

1.587899 

0.8032787 

1.476507 

0.8360656 

1.396058 

0.8655738 

1.334174 

0.9344262 

1.173275 

0.9639344 

1.086638 

1 

1 

x/L 

P/Pout 

6.52E-09 

2.714189 

0.06985247 

2.647551 

0.09936067 

2.598043 

0.1354262 

2.554725 

0.1682131 

2.517594 

0.1977213 

2.480464 

0.2337869 

2.424768 

0.2665738 

2.375261 

0.2993607 

2.33813 

0.3354262 

2.282435 

0.3682131 

2.239116 

0.401 

2.189609 

0.4370656 

2.133913 

0.4665738 

2.084406 

0.5026394 

2.03871 

0.5321475 

1.989203 

0.5682131 

1.939696 

0.6337869 

1.815928 

0.6665738 

1.754044 

0.6993607 

1.698348 

0.7354263 

1.627899 

0.7747705 

1.578391 

0.8042787 

1.516507 

0.8370656 

1.436058 

0.8665738 

1.374174 

0.8993607 

1.293725 

0.9354262 

1.213275 

0.9649344 

1.126638 

1 

1 
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Table  14.  Normalized  friction  coefficients  comparison  (see  Figure  31). 


Dh,  /zm 

C  ic 

Numerical 

C' ic 

From  eq.  (29) 

%  Difference 

1 

0.622 

0.625 

0.586 

2 

0.765 

0.769 

0.638 

10 

0.941 

0.943 

0.124 

20 

0.969 

0.971 

0.065 

40 

0.985 

0.985 

0.053 

50 

0.987 

0.988 

0.044 

60 

0.989 

0.991 

0.033 

80 

0.992 

0.993 

0.026 

100 

0.994 

0.995 

0.022 

Table  15.  Mass  flow  rates  at  various  pressure  ratios  (See  Figure  33). 


P  ir/P out 

Mass  Flow  Rate,  kg/s 

Compressible  Slip  flow 

Compressible  No-slip 

Incompressible  No-slip 

1.5 

2.93564E-13 

1.91068E-13 

1.52843E-13 

2 

6.67737E-13 

4.58537E-13 

3.05742E-13 

2.5 

1.12066E-12 

8.02481E-13 

4.58621E-13 

3 

1.65130E-12 

1.22287E-12 

6.29309E-13 

3.5 

2.58108E-12 

1.71934E-12 

7.64255E-13 

4 

2.94287E-12 

2.29230E-12 

9.17130E-13 

4.5 

3.70404E-12 

2.94126E-12 

1.06999E-12 

5 

4.54195E-12 

3.66662E-12 

1.22277E-12 

Table  16.  Average  Knudsen  number  with  varying  pressure  ratio  (See  Figure  35). 


Pin/Pout 

Average  Kn  number 

Compressible  Slip  flow 

Compressible  No-slip 

Incompressible  No-slip 

1.5 

4.45E-02 

0.044751 

0.061373 

2 

3.48E-02 

0.034909 

0.061373 

2.5 

2.85E-02 

0.028508 

0.061373 

3 

2.41E-02 

0.024046 

0.061373 

3.5 

2.09E-02 

0.020772 

0.061373 

4 

1.84E-02 

0.018271 

0.061373 

4.5 

1.64E-02 

0.016302 

0.061373 

5 

1.48E-02 

0.014712 

0.061373 
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Table  17.  Average  Mach  number  with  varying  pressure  ratio  (See  Figure  36). 


Pin/Pout 

Average  Mach  number 

Compressible  Slip  flow 

1.5 

0.002533 

0.001927 

0.002093 

2 

0.004664 

0.003665 

0.004186 

2.5 

0.006599 

0.005316 

0.006279 

3 

0.008421 

0.006917 

0.008626 

3.5 

0.010168 

0.008488 

0.010463 

4 

0.011868 

0.01004 

0.012556 

4.5 

0.013531 

0.011578 

0.014648 

5 

0.015168 

0.013107 

0.01674 

Table  18.  Cc  versus  varying  pressure  ratio,  Cc  =  fReComp/fReincomp  (See  Figure  37). 


Pin/Pout 

Q.  —  fRCComp/fR^lDcomp 

Compressible  Slip  flow 

ssnsi 

1.5 

0.693295652 

1.074244047 

1 

2 

0.757747981 

1.110755759 

l 

2.5 

0.804152495 

1.131335386 

l 

3 

0.839292896 

1.143927117 

1 

3.5 

0.863271638 

1.15194438 

1 

4 

0.889292205 

1.157532565 

1 

4.5 

0.908009958 

1.161366766 

1 

5 

0.923859309 

1.164120215 

l 

Table  19.  Comparison  of  2-D  and  3-D  results  for  Cc  -  fReComp/fReincomp  (See  Figure  44). 


Pin/Pout 

Cc  —  fRCComp/fRClncomp 

3-D  Compressible  Slip 

3-D  Incompressible  No-slip 

1.5 

0.693296 

1 

2 

0.757748 

1 

0.654161488 

2.5 

0.804152 

1 

3 

0.839293 

1 

0.759835223 

3.5 

0.863272 

1 

4 

0.889292 

1 

0.790601968 

4.5 

0.90801 

1 

5 

0.923859 

1 

0.813153222 

Experimental  Study  and  Numerical  Modeling  of  Gas  Flow  in  Microchannels  and  Micronozzles 


105 


HCET-2003-M018-001-04 


Table  20.  The  2-D  and  3-D  flow  U  -velocity  comparison  (See  Figure  41). 


Pin  /  P out 

2-D 

(U  m/s) 

3-D 

(U  m/s) 

%  Difference 

2.0 

6.55 

3.79 

42.2 

3.0 

14.89 

9.01 

39.49 

4.0 

26.50 

16.12 

39.17 

5.0 

40.81 

24.24 

40.61 

Table  21.  Nozzle  exit  thrust:  comparison  of  DSMC,  experiments,  and  Navier-Stokes  solution  (see 
Figure  54). 


RCthroat 

Exit  Thrust  (Tn),  N 

Experiment 

DSMC 

NS,  slip 
Tn=MV 

NS,  slip 
Tn=MV+PA 

0.16 

4.54E-07 

2.38E-06 

1.07E-07 

4.91E-05 

0.22 

8.82E-07 

3.21E-06 

1.83E-07 

4.92E-05 

0.38 

2.01E-06 

5.73E-06 

4.99E-07 

4.95E-05 

0.55 

2.76E-06 

7.99E-06 

9.00E-07 

4.99E-05 

0.77 

4.45E-06 

1.13E-05 

1.56E-06 

5.06E-05 

1.04 

6.45E-06 

1.53E-05 

2.62E-06 

5.17E-05 

2.64 

1.72E-05 

3.97E-05 

1.15E-05 

6.06E-05 

5.16 

3.92E-05 

7.86E-05 

3.19E-05 

8.09E-05 

7.58 

6.67E-05 

1.15E-04 

5.65E-05 

1.06E-04 

10.93 

1.13E-04 

1.66E-04 

9.52E-05 

1.44E-04 

14.45 

1.69E-04 

2.19E-04 

1.42E-04 

1.91E-04 

18.08 

2.32E-04 

2.90E-04 

1.93E-04 

2.45E-04 

22.69 

3.19E-04 

3.54E-04 

2.55E-04 

3.36E-04 

Table  22.  DSMC  simulation  results  for  the  effect  of  the  inlet  temperature  of  the  propellant  gas  at  P0  = 
1  Torr,  Dtiiroai=  1  mm  (see  Figure  56). 


R^throat 

Tn,  N 

T0  =  300  K 

T0=  1300  K 

5 

7.61E-05 

0.000226 

10 

0.000153 

0.000461 

20 

0.00031 

0.000939 

50 

0.000776 

0.002992 

100 

0.001605 

0.005945 
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Table  23.  MatLab  Image  Correlation  Velocimetry  script. 


%  ncc-normxcorr2/ image  crop/x-y  SW  size/x-y  max  distance/x-y  overlap/ 

%  parameters :  swx, swy, dxlm, dxrrn, dytm, dybm, lapx, lapy, tresh, cropimg,maporf ig 
%  19  Jan  2005 

%  Author:  Ali  K.  Reyhanogullari 
clc  ; 

%%%  CONSTANT  DEFINITIONS 
%  display  preferences 

version  =  'ICVncc23';  qscale  =  0.75;  disp_opt  =  3; 
xmmstep  =  100;  ymmstep  =50; 

%  area  to  be  processed 

tope  =  0.3;  bottomc  =  0.7;  leftc  =  0.00;  rightc  =  1.00; 

%  sample  vector  value 
sampleval  =  [20.0  0.0]; 

%%%  PRE-PROCESSING 

%  extract  raw  image  pair  and  experiment  data  from  dataset 
imgl  =  FMInput {1,2}. Img . Imgl ; 
img2  =  FMInput {1,2} . Img . Img2 ; 

%  experimental  setup  parameters 

dt  =  FMInput{l, 1} .TimeBtwPulses; 

sfactor  =  FMInput {1, 1} .CameraSetup. ScaleFactor; 

pxpitch  =  FMInput {1, 1} .CameraSetup. PxPitch; 

%  extract  parameters 

params  =  str2num(ParamStr) ; 

swx  =  params ( 1 )  ; 

swy  =  params (2); 

dxlm  =  params ( 3 ) ; 

dxrm  =  params ( 4 )  ; 

dytm  =  params  ( 5 )  ; 

dybm  =  params ( 6 )  ; 

lapx  =  params (7); 

lapy  =  params  (8); 

tresh  =  params ( 9 )  ; 

crop img  =  params  (10); 

maporfig  =  params (11); 

%  raw  image  width  and  height  in  pixels 
[n,m]  =  size(imgl); 

%  if  desired,  crop  raw  image  to  process  only  a  portion 
if  (cropimg  ==  1) 

imgl  =  imerop ( imgl , (leftc*m  topc*n  (rightc-lef tc) *m  (bottomc-topc) *n] ) ; 
img2  =  imerop (img2 , [leftc*m  topc*n  (rightc-lef tc) *m  (bottomc-topc) *n] ) ; 

end 

%  updated  image  width  and  height  in  pixels 
[nc,mc]  =  size(imgl); 

%  image  dimensions  in  mm  for  display 
xmmloc  =  [0 :xmmstep:mc] ' ; 
ymmloc  =  [0:ymmstep:nc] ' ; 

xmm  =  num2str (round (xmmloc*pxpitch(l) *sfactor*le+5) *le-2) ; 
ymm  =  num2str (round (ymmloc*pxpitch (2 ) *sfactor*le+5 ) *le-2 ) ; 

%%%  IMAGE  ENHANCEMENT 
%  contrast  scaling 
imgl  =  double ( imgl ) ; 
img2  =  double (img2) ; 

%  find  minimum  and  maximum  pixel  values _ _ 
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a  =  min ( imgl (:)),- 
b  =  max { imgl ( : ) ) ; 

%  rescale  image  matrix  according  to  these  values,  resulting  pixel  values 
%  might  be  floating  numbers  and  should  be  converted  to  integers 
if  (b  <  256) 

imgtype  =  '8-bit'; 
imgl  =  uint8 (255/ (b-a) * (imgl-a) ) ; 
img2  =  uint8(255/ (b-a) * (img2-a) ) ; 
elseif  (b  <  65536) 

imgtype  =  '16-bit'; 

imgl  =  uintl6 (65535/ (b-a) * (imgl-a) )  ; 
img2  =  uintl6 (65535/ (b-a) * (img2-a) ) ; 

else 

error ('Image  is  not  a  MATLAB  indexed  8  or  16  bit  type.'); 

end 

%%%  NCC  FOR  EACH  INTERROGATION  WINDOW  (SW) 
tic; 

%  number  of  SWs  in  x  and  y  direction  of  image  including  overlaps 
lapx  =  floor (lapx*swx) ; 
lapy  =  floor (lapy*swy) ; 

nswx  =  floor ( (me- (dxlm+dxrm) -lapx) / (swx-lapx) ) ; 

nswy  =  floor ( (nc- (dytm+dybm) -lapy) / (swy-lapy) ) ; 

sx  =  floor ( (me- (dxlm+dxrm) - (nswx*swx- (nswx-1) *lapx) ) /2 ) ; 

sy  =  floor ( (nc- (dytm+dybm) - (nswy*swy- (nswy-1) *lapy) ) /2 ) ; 

treshav  =  0; 

for  i=l:nswy 

%  upper  edge  of  current  SW 

yssw  =  sy+dytm+l+ (i-1) *swy- (i-1) *lapy; 

for  j=l: nswx 

%  left  edge  of  current  SW 

xssw  =  sx+dxlm+l+ ( j-1) *swx- ( j  — 1 ) *lapx; 

%  crop  current  SW  and  RW 

SW  =  imcrop( imgl, [xssw  yssw  (swx-1)  (swy-1)]); 

RW  =  imcrop(img2, [xssw-dxlra  yssw-dytm  (dxlm+swx-l+dxrm) . . . 

(dytm+swy-l+dybm) ] ) ; 

NCC  =  normxcorr2 (SW, RW) ; 

%  since  normxcorr2  gives  "extrapolated"  coefficients,  discard  them 
NCC  =  imcrop(NCC, [swx  swy  (dxlm+dxrm)  (dytm+dybm) ]) ; 

[NCCmax, iNCCmax]  =  max(NCC(:)); 

%  assign  shift  values,  assign  NaN  if  no  satisfactory  correlation 
if  (NCCmax  >  tresh) 

[yNCCmax, xNCCmax]  =  ind2sub(size (NCC) , iNCCmax(l) ) ; 
xshift(i,j)  =  xNCCmax- (dxlm+1) ; 
yshift(i,j)  =  yNCCmax- (dytm+1) ; 

else 

xshift(i,j)  =  0; 
yshif t ( i , j )  =  0; 

end 

%  position  for  these  vectors 
xpos(i,j)  =  xssw  -  1  +  floor (swx/2) ; 
ypos(i,j)  =  yssw  -  1  +  floor (swy/2) ; 

%  velocities  in  m/sec  calculated  using  delta_t,  pixel _pitch  and 
%  scale_factor 
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U(i,j)  =  (xshif t ( i, j ) *pxpitch (1) ) *sfactor/dt ; 

V(i,j)  =  (yshift (i, j) *pxpitch(2) ) *sfactor/dt; 

treshav  =  treshav+NCCmax; 

end 

end 

tool  =  toe; 

treshav  =  treshav/ (nswx*nswy)  ; 

Uav  =  sum(sum(U) ) / (nswx*nswy)  ; 

%%%  SEND  RESULTS  TO  FLOWMANAGER 
if  (maporfig  ==  0) 

format  compact; 

disp ( 'Analysis  with:  ');  disp (version) ;  disp('  '); 

%  flip  vector  map  vertically 

%  RECALL:  FlowManager  takes  bottom  of  the  image  as  the  vertical  0 
%  whereas  MATLAB  takes  it  as  the  top  of  the  image) 
ypos  =  nc-ypos ; 

V  =  -V; 

FMOutVector (xpos , ypos ,  U,  V)  ; 
treshav 

%  nswx,nswy,Uav,Vav, treshav, delta_t, scale_factor,pixel_pitch 
else 

%  clear  any  previous  figure 
elf  ; 

%  color  schemes 
if  (disp_opt  ==  0) 

map  =  'white';  color  =  '  b‘; 
elseif  (disp_opt  ==  1) 

map  =  'bone';  color  =  'c'; 
elseif  (disp_opt  ==  2) 

map  =  ' summer ' ;  color  =  ' k ' ; 
else 

map  =  'jet';  color  =  'r'; 

end 

cap  =  'Analysis  with  %s:  %s  Vector  Map  [%dx%d]  found  in  %4.2fs.'; 
info  =  'swx=%d  swy=%d  dxlm=%d  dxrm=%d  dytm=%d  dybm=%d'; 
info  =  streat (info, '  lapx=%d  lapy=%d  tresh=%l . 2f ’ ) ; 

info  =  streat (info, '  Uav=%2.2fm/s  treshav=%l .2f  SF=%1.4f  dt=%1.2e'); 
%  sample  vector 

U(l,nswx)  =  sampleval (1) ;  V(l,nswx)  =  sampleval (2) ; 
sample  =  '  %2.2f  m/s  \n'; 

sample  =  sprintf (sample, sqrt (sampleval (1) ^2+sampleval (2) ~2) ) ; 

hfigure  =  figured); 

set (hfigure, 'Visible ' , ' off ’ ) ; 

cap  =  sprintf (cap, version, '  ' ,nswx,nswy, tocl) ; 

info  =  sprintf ( inf o, swx, swy, dxlm, dxrm, dytm,dybm, lapx, lapy, tresh, .. . 

Uav, treshav, s factor, dt) ; 
imshow ( img2 ) ;  colormap (map ) ; 
hold  on;  grid  on;  axis  on; 

text (xmmloc-mc*0 . 015 , nc*l . 12*ones (size (xmmloc, 1) , 1) , xmm) ; 

text ( -mc*0 . 10*ones (size (ymmloc, 1) , 1) , ymmloc,ymm) ; 

text ( -mc*0 .10, -nc*0 . 06, 'mm' ) ;  text ( -mc*0 . 04, -nc*0 . 06 , 'pixel ’ ) ; 

text (me *1. 015, nc *1.12, 'mm' ) ;  text (mc*l .015,nc*1.06, 'pixel ' ) ; 

hquiver  =  quiver (xpos, ypos , U,V, qscale, color) ; 

setfhquiver, 'LineWidth ’ , 1 . 4) ; 

htext  =  text (xpos (l,nswx) ,ypos (1 ,nswx) *0 . 90 , sample) ; 
set(htext, 'EdgeColor' , 'k' , 'Color' , 'k' , 'Margin' ,nc*0.01) ; 
title (cap) ;  text (0, nc*l . 18, info) ; 

_ plot (xpos (1 ,nswx) ,ypos (1 , nswx) , 'g . ' ,  'LineWidth ' , 5)  ; _ 
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